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SUMMARY 


The  atmospheric  processes  determining  the  late-time  evolution  of  near-surface 
dust  concentrations  from  a  nuclear  burst  are  examined.  Late-time  dispersion  models, 
which  describe  the  dust  cloud  during  the  period  from  cloud  stabilization  to  many  hours 
after  initiation,  have  principally  focused  on  wind  transport  and  fallout  phenomena. 
However,  some  applications  require  knowledge  of  the  environment  near  the  ground  and 
other  processes  can  play  an  important  role  in  this  region. 

The  atmosphere  near  the  ground  is  usually  turbulent;  the  wind  blowing  over  the 
rough  surface  produces  shear-driven  turbulence,  and  daytime  heating  of  the  surface 
generates  buoyancy-driven  turbulence.  Vertical  diffusion  of  dust  particles  is  therefore  a 
significant  effect  with  vertical  velocity  fluctuations  of  the  order  of  lms~^  maintaimng 
small  particles  aloft  for  long  periods,  and  strongly  influencing  the  vertical  distribution  of 
dust.  The  turbulent  layer,  known  as  the  planetary  boundary  layer  (PBL),  varies  in  depth 
from  a  few  hundred  meters  up  to  2km  or  more  during  the  day,  but  is  usually  much 
shallower  under  nocturnal  conditions.  The  effects  of  the  PBL  are  therefore  strongly 
dependent  on  the  time  of  day  and  meteorological  conditions,  and  cannot  be  properly 
represented  by  a  fixed  diffusivity. 

The  research  reported  here  seeks  to  provide  representations  of  the  boundary  layer 
processes  for  late-time  dispersion  models.  The  parameterization  of  the  PBL  turbulence 
and  diffusion  is  discussed,  and  the  effects  of  radiative  absorption  by  the  dust  are  also 
considered.  The  deposition  of  particles  on  the  surface  through  turbulent  processes,  as 
distinct  from  gravitational  settling,  are  examined.  This  process,  known  as  dry  deposition, 
is  important  for  small  particles  with  diameters  less  than  about  30|jm,  since  the  turbulent 
deposition  rates  can  dominate  the  overall  removal  rate.  Parameterization  schemes  exist  in 
the  literature,  describing  the  effects  of  vegetative  canopies  and  rough  surfaces,  and  a 
general  representation  for  use  in  late-time  models  is  constructed.  The  structure  of  the 
PBL  is  considered  in  the  representation,  and  a  new  description  of  the  turbulent  velocity 
fluctuations  for  light  wind  conditions  is  obtained. 

The  existing  descriptions  of  the  PBL  are  generally  restricted  to  homogeneous,  i.e., 
flat,  conditions,  whereas  most  land  regions  of  the  globe  are  not  flat.  A  significant  part  of 
the  research  reported  here  was  devoted  to  extending  our  understanding  of  boundary  layer 
processes  over  complex  terrain,  and  a  large  number  of  detailed  turbulent  simulations  were 
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performed  for  a  range  of  terrain  shapes.  Results  were  obtained  for  the  surface  flux  of 
momentum,  which  determines  the  mean  boundary  layer  transport  velocity,  and  statistics 
for  the  turbulent  fluctuations  were  also  examined.  A  general  parameterization  for  the 
mean  surface  drag  and  boundary  layer  wind  is  obtained,  which  fits  numerical  results  for  a 
wide  range  of  meteorological  conditions  and  terrain  shapes.  The  representation  utilizes 
the  r.m.s.  terrain  slopes  and  elevation  variations,  and  is  constructed  with  a  proper  tensor 
form  so  that  directional  effects  are  correctly  described  for  terrain  with  dominant  ridge- 
valley  axes.  Steep  slopes  induce  lee  separation,  and  the  dispersion  of  particles  over 
terrain  was  found  to  be  strongly  enhanced  by  the  large  mean  velocity  gradients. 
Turbulence  levels  are  also  increased  by  the  shear  distortions,  and  simple  estimates  for  use 
in  late-time  dispersion  models  are  suggested.  In  contrast  to  the  dispersion  effects,  the 
turbulent  deposition  of  dust  particles  was  much  less  sensitive  to  terrain  variations.  This  is 
because  the  deposition  is  dependent  on  the  surface  shear  stress,  while  a  large  component 
of  the  force  on  the  terrain  is  due  to  the  pressure  imbalance.  The  shear  stress  is  increased 
over  hill  crests  but  is  reduced  in  sheltered  regions,  so  that  the  net  change  is  generally 
small. 
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3-4  Concentration  profiles  after  a  morning  release  into  a  convective  boundary 
layer  with  2ms“i  wind  speed  showing  the  radiative  effect  of  dust  loading. 

Solid  line  is  high  loading,  long  dashes  medium  loading,  short  dashes  low 
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effects . 

3-5  Potential  temperature  profiles  after  a  morning  release  into  a  convective 
boundary  layer  with  2ms-i  wind  speed  showing  the  effect  of  dust  loading. 

Units  are  "C,  and  the  line  patterns  are  as  in  Figure  3-4 .  25 

3-6  Vertical  velocity  variance  profiles  after  a  morning  release  into  a 
convective  boundary  layer  with  2ms"'^  wind  speed  showing  the  effect  of 
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3-7  Sensible  heat  flux  profiles  after  a  morning  release  into  a  convective 
boundary  layer  with  2ms-i  wind  speed  showing  the  effect  of  dust  loading. 

Units  are  °Kms-^  and  the  line  patterns  are  as  in  Figure  3-4 .  27 
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4-7  Profiles  of  turbulent  velocity  fluctuations  for  three  forest  types 
normaUzed  by  the  value  above  the  canopy  from  Amiro  (1990)  compared 
with  one-dimensional  canopy  model  and  LES . 
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4- 9  Probability  density  function  for  the  instantaneous  friction  velocity. 

(a)  variation  with  roughness  length  for  free  convection  conditions, 

(b)  variation  with  geostrophic  wind  speed,  dimensionless  roughness  of 

10-4  except  where  indicated . 
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5-3  Mean  flow  over  2km  ridges  with  slope  0.5  and  zo=lni.  (a)  mean 
streamlines,  contour  interval  is  200m2s-i  for  positive  values,  lOm^s  ^  for 
negative  values  (shown  dashed),  (b)  mean  transverse  velocity 
component,  contour  interval  of  . . 
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5-5  Time  history  of  the  domain-averaged  surface  forces  for  flow  over  2km 
ridges  with  slope  0.5  and  zo=lm.  Forces  are  displayed  per  unit  horizontal 
surface  area  and  per  unit  fluid  density . 
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SECTION  1 


INTRODUCTION 


1.1  BACKGROUND. 

The  fate  of  dust  particles  injected  into  the  atmosphere  by  nuclear  explosion  has 
been  studied  extensively.  Most  of  the  studies  have  concentrated  on  high  altitude  effects, 
where  small  particles  can  persist  for  extended  periods  with  possible  climatic  impact. 
Long  range  radioactive  fallout  calculations  are  also  principally  concerned  with 
atmospheric  dispersion  phenomena  in  the  upper  atmosphere.  However,  there  are  a 
number  of  applications  where  the  evolution  of  the  dust  cloud  near  the  ground  is  of 
interest.  For  example,  aircraft  missions  include  extensive  flight  periods  at  low  altitude, 
exposing  the  aircraft  to  persistent  degradation  in  the  near-surface  environment.  The 
cumulative  effects  of  exposure  to  relatively  small  concentrations  of  dust  are  particularly 
significant,  since  this  can  cause  engine  failure  or  impair  visibility  through  transparent 
surfaces.  The  possibility  of  accumulated  damage  over  long  flightpaths  demands  the  study 
of  the  late-time  evolution  of  the  dust  field,  since  smaller  particles  can  remain  aloft  for 
many  hours. 

Large  masses  of  dust  can  be  lofted  into  the  atmosphere  by  a  nuclear  blast,  and 
there  are  two  major  paths  for  the  dust  to  appear  at  low  levels.  First,  particles  in  the  main 
cloud,  which  can  rise  many  kilometers,  will  fall  out  under  gravitational  acceleration  and 
pass  through  low  levels  at  some  later  time.  This  is  the  pathway  for  much  of  the 
radioactive  fallout  and  has  been  studied  previously.  However,  a  more  persistent  threat  is 
posed  by  the  high  density  pedestal,  formed  by  dust  which  is  lifted  from  the  surface  but  is 
not  carried  upward  with  the  main  cloud.  The  fate  of  a  pedestal  particle  depends  on  many 
factors,  both  its  own  size  and  density,  and  also  meteorological  and  terrain  factors. 

The  region  near  the  ground  surface  is  often  turbulent,  forming  the  planetary 
boundary  layer,  which  usually  extends  on  the  order  of  1000m  in  the  vertical.  The  flow  in 
the  boundary  layer  is  influenced  by  surface  processes,  including  heat  and  moisture  fluxes 
and  turbulent  drag  on  the  surface.  These  processes  induce  turbulent  motions  with  the 
boundary  layer,  and  depend  on  the  surface  conditions  as  well  as  the  solar  and  infra-red 
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radiative  fluxes.  The  turbulent  mixing  in  the  planetary  boundary  layer  can  mamtain  small 
particles  in  suspension,  while  surface  features  can  augment  the  gravitational  deposition 
process  by  absorbing  particles.  In  addition,  cloud  precipitation  processes  can  scavenge 
the  dust  from  the  atmosphere.  All  these  processes  depend  on  the  size  and  type  of  dust 
particle,  and  therefore  the  initial  pedestal  formation  is  critical,  since  this  will  depend  on 
soil  or  surface  type  as  well  as  the  burst  type.  It  is  clear  that  any  assessment  of  low-level 
dust  hazard  requires  a  reliable  description  of  both  the  pedestal  generation  mechanism  and 
the  late-time  boundary  layer  dispersion  and  deposition  mechanisms. 


Current  late-time  cloud  models  have  not  devoted  much  effort  to  the  near-surface 
dust.  The  objective  of  the  research  reported  here  is  the  development  of  algorithms  to 
allow  the  late-time  cloud  models  to  describe  the  late-time  processes  of  turbulent  diffusion 
and  deposition  in  the  planetary  boundary  kyer. 


1.2  PHYSICAL  MECHANISMS. 

Much  attention  has  been  paid  to  the  lofted  dust  cloud  from  a  nuclear  burst,  which 
can  carry  material  up  to  the  stratosphere  where  it  will  persist  for  a  very  long  time.  In 
addition  to  the  main  cloud,  however,  there  can  also  be  a  significant  'pedestal  of  dust 
swept  up  from  the  surface  by  the  initial  blast  wave,  and  also  by  the  inflow  feedmg  the 
rising  fireball.  There  is  a  potentially  large  mass  of  dust  lofted  only  into  a  shallow 
(<  100m  deep)  layer  near  the  surface,  which  does  not  become  entrained  into  the  rising 
cloud,  and  therefore  remains  close  to  the  surface  when  only  burst-induced  dynamics  are 
considered.  FoUowing  decay  of  the  burst  flow  field,  however,  the  ambient  atmosphenc 
flow  effects  will  begin  to  dominate  and  control  the  low-level  dust  evolution.  The 
following  sections  discuss  the  various  phenomena  occurring  in  the  boundary  layer  which 
affect  the  concentration  of  near-surface  dust 

1.2.1  Turbulent  Diffusion. 

The  atmosphere  near  the  surface  of  the  Earth  is  usually  in  turbulent  motion.  The 
wind  blowing  over  the  rough  surface,  and  the  buoyancy  effects  of  daytime  surface 
heating  produce  turbulent  vertical  motions  that  mix  airborne  material  both  vertically  and 
horizontally.  The  surface-generated  turbulent  layer  is  often  capped  by  a  stable 
temperature  inversion,  which  marks  the  vertical  extent  of  the  mixing  and  is  formed  as  the 
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torbulent  boundary  layer  grows  upward  into  the  stably-stratifed  atmosphere.  The  depth 
of  the  mining  region  can  vary  with  meteorological  conditions  and  time  of  day,  but  usually 
reaches  between  about  5(X)m  and  2000m  during  the  daytime.  Shallower  layers  are 
associated  with  reduced  heating,  due  to  the  presence  of  cloud,  for  example,  or  due  to 
large  scale  subsidence  overlying  air.  Under  nocturnal  cooling  conditions,  the  surface 
layer  becomes  stably-stratified  and  the  mixing  depth  and  turbulence  intensities  are 
reduced. 

Vertical  mixing  can  play  an  important  role  in  the  late-time  dispersion  of  dust 
clouds,  as  demonstrated  by  Sykes,  Parker  and  Henn  (1993).  Daytime  turbulent  velocity 
fluctuations  are  on  the  order  of  Ims-i,  and  can  therefore  dominate  the  vertical  transport  of 
particles  with  gravitational  settling  speeds  smaller  than  this.  Thus,  particles  smaller  than 
about  SOpjn  will  be  strongly  affected  by  the  mixing  and  can  remain  aloft  for  much  longer 
than  under  quiescent  conditions.  Many  schemes  have  been  proposed  to  represent  the 
diffusion  processes  in  the  atmospheric  boundary  layer,  and  appropriate  methods  for  late¬ 
time  dispersion  modeling  will  be  discussed  in  Section  2,  and  the  radiative  effects  of  a 
dust  cloud  on  the  turbulent  boundary  layer  are  considered  in  Section  3.. 

1.2,2  Dry  Deposition. 

The  absorption  of  a  contaminant  substance  at  the  surface  is  known  as  dry 
deposition.  Dry  deposition  can  include  moisture  effects,  but  is  distinguished  from  the 
direct  deposition  of  dust  in  precipitation  such  as  rain  or  snow.  The  process  is  complex 
and  is  controlled  by  a  number  of  mechanisms.  Idealized  relationships  have  been 
established  for  massless  scalar  contaminants  in  terms  of  the  surface  conditions  but  these 
relationships  are  not  directly  applicable  to  the  dust  particle  deposition.  The  finite-size 
particles  do  not  precisely  follow  the  air  flow  and  they  do  not  diffuse  with  the  same 
molecular  diffusivity  as  a  scalar  contaminant  such  as  a  trace  gas.  The  deposition  is 
usually  modeled  in  terms  of  a  sum  of  resistances  for  each  of  the  various  processes,  such 
as  turbulent  transfer  or  molecular  diffusion;  a  rapid  transfer  rate  is  represented  by  a  low 
resistance.  Two  general  types  of  surface  are  considered,  a  smooth  surface  and  a 
vegetative  canopy. 

In  the  atmospheric  context,  a  smooth  surface  is  represented  by  surfaces  such  as 
desert  sand  or  snow.  The  surface  conditions  are  obtained  from  the  assumption  that 
particles  which  impact  the  surface  will  stick.  The  flux  of  material  on  to  the  surface 
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therefore  depends  on  the  rate  of  impact.  Very  small  particles  tend  to  follow  the  air 
around  surface  features,  while  large  particles  simply  fall  rapidly  to  the  ground. 
Intermediate  sizes  wUl  be  unable  to  follow  the  air  around  small  surface  features  and  will 
be  deposited  on  those  features.  Small  particles  can  also  reach  the  surface  by  Brownian 
motion  if  they  are  brought  sufficiently  close.  The  deposition  is  therefore  governed  by  the 
particle  response  time,  the  turbulence  acceleration  rates,  the  scale  of  the  surface  features, 
and  the  Brownian  diffusion  rate. 

A  vegetative  canopy  is  distinguished  by  its  three-dimensional  nature,  i.e.,  there  is 
a  significant  flow  through  the  canopy.  A  surface  covered  by  trees  in  full  leaf  can  provide 
a  large  particulate  flux  onto  the  surface,  since  the  dust  in  the  canopy  layer  will  be  swept 
past  a  vast  number  of  individual  leaves  greatly  increasing  the  chance  of  impact.  The 
actual  deposition  process  is  extremely  complex,  reflecting  the  microstructure  of  the  leaf 
structure  of  the  various  species  in  the  canopy.  The  flux  also  depends  on  the  rate  of  which 
the  air  sweeps  past  the  leaves,  and  this  is  controlled  by  the  drag  effects  of  the  canopy.  It 
is  well-known  that  a  dense  canopy  reduces  the  wind  speed,  so  all  these  effects  need  to  be 
modeled  if  we  are  to  determine  the  fate  of  low  level  dust. 

Representations  of  the  deposition  process  will  be  discussed  in  Section  4. 

1.2.3  Wet  Deposition  and  Scavenging. 

The  condensation  of  water  vapor  into  liquid  droplets  or  frozen  ice  particles 
introduces  a  new  range  of  mechanisms  for  dust  removal.  The  basic  phenomenon  requires 
the  incorporation  of  the  dust  particle  into  a  water/ice  droplet,  which  can  then  fail  to  the 
ground  and  remove  the  dust  from  the  atmosphere.  The  mechanisms  for  coalescence  are 
manifold,  including  collision  impact,  nucleation,  molecular  diffusion,  electrical  effects,  or 
phoretic  effects  due  to  local  gradients  of  molecular  properties. 


The  precipitation  scavenging  process  involves  hydrometeors  falling  through  the 
dusty  boundary  layer,  collecting  dust  particles  and  depositing  them  on  the  surface.  In 
general,  the  scavenging  mechanisms  is  very  efficient  in  removing  small  particles,  so  that 
a  heavy  rain  can  'clean'  the  boundary  layer  relatively  quickly  (Schwartz,  1992).  The 
collection  efficiency  does  depend  on  the  relative  sizes  of  the  dust  and  the  precipitation, 
and  is  different  for  different  types  of  precipitation,  e.g.  rain,  snow,  hail. 
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In  general,  we  do  not  expect  the  meteorological  input  to  late-time  dispersion 
models  to  contain  information  on  precipitation  rates.  Precipitation  is  difficult  to  predict 
and  is  usually  very  intermittent.  It  is  therefore  inappropriate  to  attempt  a  detailed 
representation  of  the  microphysical  processes  involved  in  precipitation  scavenging,  and  a 
much  simpler  approach  is  recommended.  If  the  meteorological  input  contains 
information  on  the  precipitation  fields,  then  the  precipitation  fall  speed  can  simply  be 
added  to  the  particle  faU  speed.  This  assumes  a  perfectly  efficient  collection  mechanism, 
i..e.,  every  dust  particle  is  immediately  attached  to  a  droplet,  but  there  is  evidence  that 
small  particles  are  indeed  collected  efficiently,  and  large  particles  will  already  have  a 
significant  fall  speed.  At  present,  this  simple  treatment  of  precipitation  scavenging  is  all 
that  can  be  justified  in  light  of  the  available  meteorological  information. 

1.2.4  Terrain  Effects. 


The  earth's  surface  is  generally  non-uniform,  with  changes  in  surface  type  and 
elevation  on  all  scales.  The  processes  described  in  the  previous  sections  are  all  local 
phenomena,  and  will  respond  to  local  features  such  as  hills  or  forests.  Changes  in  surface 
roughness,  e.g.,  wheat  field  to  forest,  are  generally  significant  only  near  the  surface  and 
can  be  modeled  using  purely  local  analysis,  although  the  non-equilibrium  situations  in 
regions  of  such  change  can  result  in  large  effective  deposition  rates.  For  example,  the 
wind  speed  takes  some  time  to  slow  down  after  entering  a  forest  canopy  and  there  is 
greatly  enhanced  deposition  in  the  forest  edges  where  the  winds  are  higher  than 
equilibrium. 


In  contrast  to  land  use  variations,  topographic  features  induce  large  changes  in  the 
entire  boundary  layer  and  therefore  require  more  careful  study.  One  of  the  main  effects 
of  topography  is  to  cause  the  near-surface  wind  to  speed  up  over  the  hiU  tops.  This  can 
increase  dry  deposition  in  these  regions,  and  this  may  not  be  balanced  by  a  corresponding 
decrease  in  the  valley  regions.  The  hills  also  affect  the  boundary  layer  turbulence, 
modifying  the  rate  at  which  dust  is  diffused  vertically. 

The  representation  of  complex  terrain  effects  is  discussed  in  Section  5. 
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SECTION  2 

PLANETARY  BOUNDARY  LAYER 

One  of  the  most  important  mechanisms  operating  near  the  surface  of  the  earth  is 
turbulent  diffusion.  Mechanical  generation  of  turbulence  by  surface  roughness  elements 
such  as  forests,  and  buoyant  generation  from  surface  heating  during  daytime  conditions 
over  land,  produces  a  turbulent  boundary  layer  that  mixes  poUutants  throughout  its  depth. 
This  planetary  or  atmospheric  boundary  layer,  often  referred  to  as  the  PBL,  has  been 
studied  extensively  from  the  point  of  view  of  pollutant  dispersion,  and  many  of  its 
properties  are  now  weU  known.  There  are  many  references  covering  the  structure  and 
dynamics  of  the  planetary  boundary  layer,  and  its  dispersion  properties.  Detailed 
descriptions  can  be  found  in  Lumley  and  Panofsky  (1964),  Csanady  (1973),  PasquiU 
(1974),  Panofsky  and  Dutton  (1984),  and  Wyngaard  and  Venkatram  (1988). 

The  atmospheric  boundary  layer  over  land  typically  undergoes  a  strong  diurnal 
variation,  with  nocturnal  cooling  of  the  surface  suppressing  turbulent  activity  at  night. 
Surface  heating  from  solar  radiation  warms  the  lower  layers  of  the  atmosphere  in  the 
morning  and  the  turbulent  mixing  layer  deepens.  The  growth  is  rapid  durmg  the  mormng 
transition  hours  but  slower  through  the  rest  of  the  day.  The  mixing  depth  typically  ranges 
from  5(X)m  to  2000m,  depending  on  local  conditions,  and  turbulent  eddy  velocities  are 
often  a  few  meters  per  second.  The  inversion-capped  situation,  where  the  mixed  layer 
depth  is  limited  by  an  overlying  stably-stratified  air  mass,  is  very  common.  The  weather 
under  such  conditions  can  be  clear  skies,  shallow  cumulus  clouds,  or  a  low-level  stratus 
deck.  Deep  convection  conditions  represent  a  breakdown  of  the  surface  boundary  layer, 
with  low-level  air  being  carried  high  into  the  troposphere  by  the  convective  cells.  In 
some  sense,  one  can  consider  the  entire  troposphere  to  be  part  of  the  boundary  layer 

under  these  conditions. 

The  turbulence  conditions  in  the  PBL  are  classified  by  the  relative  importance  of 
buoyancy  generation,  i.e.  stable,  neutral ,  or  unstable.  In  unstable  conditions,  the  heating 
of  the  surface  is  the  dominant  process,  causing  surface  air  to  rise  in  turbulent  convection 
eddies  through  the  boundary  layer,  whUe  stable  conditions  occur  during  nocturnal  cooling 
conditions  and  the  mixing  process  expends  energy  lifting  the  heavier  surface  air.  Neutral 
conditions  are  shear-dominated,  i.e.  high  wind  speeds  or  very  smaU  buoyancy  effects,  and 
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turbulence  profiles  are  similar  to  those  found  in  wind-tunnel  boundary  layers.  The 
buoyancy  effects  are  characterized  by  the  Monin-Obukov  length,  L,  and  the  mixing  or 
boundary  layer  depth,  zi.  The  Monin-Obukov  length  is  the  height  at  which  the  buoyancy 
generation  of  turbulence  equals  mechanical  generation,  and  therefore  provides  the 
appropriate  scale  for  distinguishing  between  the  two  mechanisms.  A  similarity  theory  for 
the  surface  layer  has  been  constructed  in  terms  of  z/L,  providing  a  description  of  the 
profiles  of  various  mean  quantities  The  atmosphere  introduces  other  complicating 
factors,  such  as  the  Coriolis  effects  due  to  the  Earth's  rotation  which  produce  a  turning  of 
the  wind  direction  with  height  known  as  the  Ekman  profile.  Also  moisture  transport  and 
condensation  effects  can  release  energy  in  the  cloud  layer,  and  modify  the  turbulence 
profiles. 


Our  main  purpose  in  determining  the  structure  of  the  PEL  is  the  characterization 
of  the  mean  boundary  layer  transport  and  the  turbulent  diffusion  rates.  This  is 
accomplished  by  means  of  parameterization  schemes,  i.e.,  the  boundary  layer  is  defined 
by  a  small  number  of  parameters,  from  which  the  mean  flow  and  diffusion  rates  can  be 
deduced.  The  most  important  parameters  for  the  PEL  are 

a)  the  geostrophic,  or  free  stream  velocity,  Uq 

b)  the  surface  roughness  height,  zo 

c)  the  surface  heat  flux.  Ho 

d)  the  boundary  layer  depth,  Zi 

These  parameters  can  be  used  to  derive  two  important  velocity  scales  that  define  the 
turbulence  levels  in  the  PEL.  The  first  is  the  surface  friction  velocity,  «*,  which  is 
defined  by  the  relation 

where  T,  is  the  average  surface  stress,  and  Pa  is  the  air  density.  The  surface  stress 
depends  on  the  atmospheric  stability  as  well  as  the  surface  roughness  and  free-stream 
wind  speed.  We  shall  discuss  the  relationship  further  below. 


The  second  velocity  scale  characterizes  the  buoyancy  effects,  and  is  known  as  the 
convective  velocity  scale  (Deardorff,  1970),  w».  The  definition  is 


w* 


V^o 


,1/3 


(2.2) 
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(2.3) 


and  is  related  to  the  Monin-Obukov  length. 


through  the  equation 

£i.  ^  ^  (2.4) 

L 

Here,  k  is  von  Karman's  constant,  which  is  usually  assumed  to  be  0.4,  g  is  the 
gravitational  acceleration,  and  To  is  a  reference  temperature  for  the  Boussinesq 
approximation.  Hq  is  actually  specified  as  a  temperature  flux  in  the  above  equations,  and 
the  relative  density  flux  is  then  obtained  after  division  by  To,  which  should  represent  an 
average  temperature  for  the  boundary  layer. 

Before  discussing  the  determination  of  «♦  and  w*,  we  consider  the  representation 
of  the  turbulent  diffusion.  Our  approach  to  this  problem  has  been  based  on  second-order 
closure  techniques,  which  provide  an  estimate  of  the  diffusion  rates  in  terms  of  velocity 
fluctuation  correlations.  This  is  a  more  general  and  systematic  approach  than  empirical 
schemes  based  on  curve-fitting  experimental  dispersion  results,  but  requires  the 
specification  of  the  velocity  fluctuation  statistics. 

The  two  simplified  situations,  where  profiles  are  well  understood,  are  provided  by 
neutral  conditions  and  free  convection  conditions.  In  neutral  flow,  the  surface  heat  flux  is 
zero,  or  at  least  is  negligible  in  comparison  with  the  mechanical  shear  effects.  The 
turbulence  is  characterized  by  u*  only  in  this  situation.  Free  convection  occurs  under 
very  light  wind  conditions  with  a  positive  surface  heat  flux,  and  the  turbulence  is 
therefore  dependent  on  w*  only. 

We  represent  the  turbulence  velocity  correlations  as 

for  neutral  conditions,  and 

for  free  convection.  Here,  the  overbar  denotes  an  average  value  and  the  prime  denotes  a 
fluctuation  from  the  average. 


(2.5) 

(2.6) 
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Neutral  profiles  are  relatively  simple  linear  functions.  Wind  tunnel  studies  of  the 
aerodynamic  boundary  layer  (Klebanoff,  1955;  Townsend,  1976),  and  numerical 
calculations  of  the  neutral  Ekman  layer  (Spalart,  1989,  Mason  and  Thompson,  1987)  both 

support  velocity  variance  profiles  proportional  to  ^1  —  ^  Appropriate  dimensionless 
profiles  are 


^..=5(1-^,) 

2.5(1 

f33=l-S(l-X) 


Free  convection  profiles  are  available  from  the  laboratory  experiments  of  Deardorff 
(1970),  and  LES  calculations  (Mason,  1989;  Schmidt  and  Schumann,  1989),  and  a 
reasonably  good  representation  is  provided  by  the  expressions 


Gii  =  G22=  0.13(1  + 1.5^^"') 


These  representations  approximate  the  turbulence  profiles  under  the  idealized  conditions, 
and  the  value  for  general  conditions  can  be  estimated  as  the  sum  of  the  two  components. 
For  stable  conditions,  the  convective  velocity  scale  is  zero  and  the  neutral  profile  shape  is 
used.  However,  the  surface  friction  velocity  will  generally  be  reduced  by  the  presence  of 
stable  stratification  near  the  surface,  and  the  boundary  layer  depth,  z,-,  is  estimated  as  5L. 


The  vertical  diffusivity,  Kf,  can  be  estimated  from  the  turbulence  profiles  using 
the  equilibrium  version  of  the  second-order  closure  equations  (Lewellen,  1977).  This 
gives 


K,= 


To  q  bs 


Aq 


1  g  Ar  dd 


1  + 


lAbs  Tq  (f  dz 


(2.9) 
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where  and  9  is  the  potential  temperature.  The  heat  flux  profile  can  be 

represented  as  a  linear  function 

VF  =  Ho(i- 

and  the  turbulence  length  scale  is  modeled  very  simply  as 

1  ^  1  I  i _  (2.11) 

a2  (0.3zif  {0.65zf 

providing  a  transition  from  the  linear  behavior  near  the  surface  to  a  constant  value  m  the 
mixed  layer.  The  turbulence  model  constants,  A,  b,  and  s,  take  the  values  0.75,  0.125, 
and  1.8,  respectively. 


The  potential  temperature  gradient,  can  be  assumed  to  be  zero  for  neutral  and 

oz 

convective  conditions,  since  the  temperature  is  well-mixed  except  in  a  thin  superadiabatic 
layer  at  the  surface.  For  stable  conditions,  the  potential  temperature  gradient  must  be 

estimated  separately. 

Several  examples  of  the  vertical  diffusivity  profiles  from  (2.9)  are  shown  for 
different  values  of  and  w*  in  Figure  2-1.  The  profiles  are  made  dimensional  by 
assuming  an  inversion  depth  of  1000m.  The  profile  shape  shows  an  elevated  maximum, 
since  Kt  goes  to  zero  at  the  ground,  where  the  length  scale  vanishes,  and  at  the  inversion, 
where  the  turbulence  levels  are  low.  For  free  convection  conditions,  m*=0,  the  diffusivity 
profile  is  in  good  agreement  with  the  values  obtained  by  Wyngaard  and  Brost  (1984) 
from  Large-Eddy  Simulations;  the  simplified  representation  (2.9)  is  not  able  to  reproduce 
the  asymmetry  between  upward  and  downward  diffusion,  but  the  predicted  magnitude  is 
representative  of  the  average  value.  The  profiles  show  an  interesting  variation  as  the 
wind  shear  is  increased  from  the  free  convection  situation  with  no  wind.  The  diffusivity 
is  initially  reduced  as  u*  increases,  in  spite  of  the  fact  that  the  turbulence  intensities  are 
additive  in  (2.9)  and  therefore  must  monotonically  increase  with  u*  .  This  anomalous 
behavior  in  the  diffusivity  is  actually  observed  in  the  Large-Eddy  Simulations  of  Mason 
(1992),  where  it  is  attributed  to  the  reduction  in  turbulence  scale  in  the  presence  of  wind 
shear.  We  do  not  modify  the  definition  of  A  in  our  simple  parameterization,  but  the 
reduction  in  Kt  is  due  to  a  reduction  in  the  turbulence  timescale,  A/q,  which  reduces  the 
contribution  from  the  heat  flux  term  in  (2.9).  It  is  encouraging  that  the  closure  model 
prediction  reproduces  this  observed  shear  phenomenon. 


10 


U' 


0 


100  200 
Diffusivity,  iTj’,  (m^s'i) 


300 


Diffusivity, /sTj-,  (m^s 


Figure  2-1.  Vertical  diffusivity  profiles  for  several  combinations  of  surface  friction 
velocity  and  convective  velocity  from  (2.9).  (a)  w*=lms“^,  solid  line 
u*=0,  long  dash  M»=0.5ms“i,  short  dash  M*=lms"i;  (b)  M*=0.5ms“^  solid 
line  w*=0,  long  dash  H'»=0.5ms“^  short  dash  w»=lms“^. 


The  turbulence  and  diffusion  parameterizations  discussed  above  are  predicated  on 
the  ability  to  specify  the  boundary  layer  parameters,  u*  and  w*.  In  many  cases,  these 
parameters  will  not  be  available  from  the  meteorological  input  data,  which  might  consist 
of  a  number  of  wind  observations  or  a  gridded  three-dimensional  wind  field.  Many 
dispersion  modelers  have  faced  the  difficulty  of  characterizing  the  PEL  using  routine 
observational  data,  and  there  are  therefore  several  schemes  in  existence  to  perform  this 
task.  For  example,  the  CRDEC  dispersion  model  NUSSE  (Saucier,  1987)  contains 
algorithms  to  determine  the  PEL  turbulence  levels.  However,  we  have  previously  used 
the  more  general  scheme  incorporated  into  METPRO  (Paine,  1987),  which  is  the 
meteorological  preprocessor  developed  for  the  U.S.  Environmental  Protection  Agency. 
METPRO  is  an  extension  of  PREPRO  (Paine  and  Hanna,  1986),  which  was  developed 
for  the  Electric  Power  Research  Institute  as  part  of  the  Plume  Model  Validation  and 

Development  program. 

The  details  of  the  mathematical  modeling  in  METPRO  can  be  found  in  the  cited 
reference,  but  a  brief  outhne  of  the  methodology  is  provided  here.  The  two  surface 
parameters,  characterizing  the  momentum  and  heat  fluxes,  are  obtained  from  a  surface 
layer  analysis  using  the  Monin  Obukov  simUarity  theory.  The  surface  layer  descnption 
requires  a  surface  roughness  value.  Given  the  wind  speed  at  a  reference  height,  the 
friction  velocity,  u.,  can  be  obtained  if  the  Monin  Obukov  length,  L,  defined  in  (2.3),  is 
known.  This  requires  knowledge  of  the  surface  sensible  heat  flux,  which  is  derived  from 
a  surface  energy  balance,  if  no  micrometeorological  thermodynamic  information  is 
provided.  If  detailed  temperature  profile  data  is  available,  then  surface  layer  similarity 
theory  can  be  used  to  estimate  the  flux 

The  surface  energy  balance  is  determined  by  the  incident  solar  radiative  flux,  the 
soil  heat  flux,  and  the  turbulent  heat  and  moisture  fluxes  into  the  air.  The  surface 
moisture  flux  is  an  important  practical  consideration  since  a  large  fraction  of  the  incident 
radiative  heating  can  be  used  to  provide  the  latent  heating  necessary  to  evaporate  water 
from  the  surface.  This  is  true  not  only  over  open  water  surfaces,  but  also  over  vegetative 
canopies.  The  ratio  of  the  sensible  heat  flux  to  the  latent  heat  flux  in  the  moisture  is 
known  as  the  Eowen  ratio.  METPRO  estimates  the  incident  radiative  flux  from  the  solar 
elevation  angle,  which  is  a  function  of  latitude,  time  of  day  and  time  of  year.  The  flux 
reaching  the  ground  is  reduced  by  a  factor  dependent  on  the  cloud  cover,  if  this  is 
available  from  the  observations.  The  relative  magnitude  of  the  soil  flux  and  turbulent 
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heat  and  moisture  fluxes  is  based  on  a  parameterization  of  the  soil/canopy  moisture 
content. 
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SECTION  3 

RADIATIVE  EFFECTS 


3.1  RADIATIVE  TRANSFER  MODELS. 


As  discussed  in  the  previous  section,  the  structure  of  the  PBL  is  strongly  affected 
by  the  surface  heat  flux,  which  is  largely  determined  by  the  incoming  solar  radiation 
during  the  daytime.  It  is  clear  that  the  radiation  at  the  surface  depends  on  natural  cloud 
cover,  so  we  must  consider  the  possibility  that  a  dust  cloud  might  also  modify  the  local 
boundary  layer  through  changes  in  the  solar  radiation  reaching  the  ground. 

As  a  preliminary  assessment  of  the  radiative  transfer  properties  of  a  dust-laden 
layer,  we  consider  an  idealized  homogeneous  situation.  Suppose  we  have  a  Uon  deep 
layer  of  monodisperse  dust  particles  of  diameter,  dy  and  mass  concentration  of  c.  The 
dusty  atmosphere  is  both  an  absorbing  and  a  scattering  medium  for  the  solar  radiation, 
and  we  wish  to  estimate  the  relative  fractions  of  the  solar  energy  that  are  reflected  back  to 
the  sky,/r,  absorbed  by  the  dust, /a,  or  transmitted  through  to  the  surface,/,.  The  incident 
radiative  flux  is  Fq  at  an  angle  do  to  the  vertical,  and  the  surface  albedo  is  Os,  as 
illustrated  schematically  in  Figure  3-1.  We  assume  that  the  clear  air  is  completely 
transparent  for  the  purposes  of  a  simple  calculation.  The  latter  assumption  makes  the 
vertical  location  of  the  dust  layer  irrelevant,  since  radiation  is  transmitted  without  loss 
above  and  below  the  layer.  Conservation  of  energy  requires  that 

/a  +  /,(l-«,)  +  /r=l 

i.e.,  the  total  input  from  the  sun  is  either  absorbed  by  the  dust  or  the  ground,  or  is 
reflected  back  upward.  All  the  quantities  discussed  in  this  section  refer  to  broadband 
averages  over  the  entire  spectral  region  covering  the  solar  flux.  A  detailed  treatment  of 
the  spectral  variabUity  of  the  particle  properties  is  beyond  the  scope  of  this  study,  which 
should  be  regarded  as  a  bulk  parameterization  suitable  for  practical  calculation. 

Individual  particle  properties  are  characterized  by  a  single  scattering  albedo,  (Do, 
and  a  scattering  phase  function,  where  ^  and  li'  represent  the  incident  and 

scattering  directions,  respectively.  The  albedo  specifies  the  fraction  of  incident  energy 
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Figure  3-1.  Schematic  representation  of  solar  radiative  flux  transfer  through  dusty  layer. 


scattered  by  an  individual  particle,  while  the  phase  function  gives  the  directional 
dependence  of  the  scattered  wave  function.  For  a  plane  incident  wave  encountering  a 
horizontally  homogeneous  scattering  medium,  the  radiative  intensity,  /,  satisfies  the 
transfer  equation 

=  -7  +  ^  J  dn' (3.2) 
dx  2 


where  T  is  the  optical  depth,  and  the  integral  extends  over  all  scattering  directions,  with  fj. 
and  equal  to  the  cosine  of  the  angle  between  the  wave  direction  and  the  vertical.  The 
optical  depth  of  the  layer  is  defined  as 


(3.3) 
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(Friediander,  1977)  where  d  is  the  particle  diameter,  N  is  the  number  density,  and  Kext  is 
the  particle  extinction  efficiency  (sometimes  called  the  attenuation  coefficient  or 
turbidity).  The  extinction  efficiency  includes  both  scattering  and  absorption  effects  of  the 
particles  and  depends  on  the  radiation  wavelength,  but  for  short  waves  relative  to  the 
particle  diameter,  Mie  theory  gives  Kext  —  2-  This  assumption  is  valid  for  solar  radiation, 
which  peaks  around  0.5p.m  wavelength,  interacting  with  particles  larger  than  1pm  and 
will  be  used  for  aU  our  subsequent  calculations.  The  quantity  t  is  sometimes  referred  to 
as  the  geometric  or  extinction  optical  depth,  and  should  not  be  interpreted  as  a  measure  of 
the  absorption  in  the  layer.  A  proper  calculation  of  the  radiative  transfer  requires 
treatment  of  the  scattering  effects  and  wiU  be  discussed  below.  The  constant  efficiency 
gives  a  simple  relation  between  the  optical  depth  and  the  mass  concentration  of  dust 
particles.  If  we  have  a  mass  loading,  c  per  unit  volume,  over  a  layer  of  depth  H,  then  the 
optical  depth  of  the  layer  is 

r  =  —  (3.4) 

Pd^ 

where  pd  is  the  particle  material  density. 

The  full  radiative  transfer  problem  is  still  extremely  complex,  even  neglecting  the 
wavelength  dependence,  since  a  real  dust  layer  involves  multiple  scattering  as  radiation 
encounters  many  particles  in  passing  through  the  layer.  Each  encounter  involves  the 
phase  function,  and  a  proper  treatment  requires  a  sophisticated  Monte-Carlo  calculation 
of  the  multiple  interaction  or  an  equivalent  technique.  We  seek  a  simplified  treatment 
that  provides  a  direct  estimate  of  radiative  transfer  for  the  multiple  scattering  regime. 

From  a  survey  of  the  atmospheric  radiation  literature,  the  delta-Eddington  two- 
stream  model  was  selected  as  the  most  appropriate  for  this  application  (Joseph  et  al., 
1976).  This  model  is  an  extension  of  the  Eddington  model  of  Shettle  and  Weinman 
(1970)  to  account  for  strong  directional  asymmetry  that  is  typically  found  in  particulate 
scattering.  Two-stream  models  represent  the  net  vertical  flux  of  broadband  radiation  as  a 
sum  of  upward,  and  downward,  fluxes.  This  approximation  is  suitable  for  quasi- 
homogeneous  situations  in  the  horizontal  direction.  The  delta-Eddington  model 
approximates  the  scattering  phase  function  by  a  delta  function  in  the  forward  direction, 
plus  a  single  spherical  harmonic,  i.e., 

P(n,ll')  =  +  (3.5) 
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where  /  is  the  fractional  scattering  into  the  forward  direction,  and  g'  is  the  asymmetry 
factor  of  the  continuous  part  of  P.  These  two  quantities  can  be  defined  in  terms  of  the 
overall  asymmetry  factor,  g,  giving 

/  =  /  (3.6a) 


8'  = 


1  +  ^ 


(3.6b) 


The  Eddington  approximation  assumes  that  the  radiance  can  be  represented  as 

/(T,/t)  =  /o(T)  +  /i/i(T) 


(3.7) 


i.e.,  an  expansion  of  the  directional  dependence  in  spherical  harmonic  functions  with 
truncation  at  first  order.  These  simplifying  assumptions  for  the  phase  function  and  the 
radiance  allow  an  analytic  integration  over  all  forward  and  backward  scattered  directions, 
as  required  for  the  two-stream  approximation,  to  give 

F'^  =  ;E(^/o(T)-|/i/i(T)j,  =  n{^Q{x)+^txIx{x'^  (3.8) 


The  expressions  for  /q  and  I\  within  the  layer  are  from  Shettle  and  Weinman 

(1970), 

/o(T)  =  (3.9a) 

h{x)  =  p[Cxe-’^^' (3.9b) 

where 

p  =  [3(l-m')/(l-5'©')f 

a  =  3(o%pI{i  +  g'{l  -  m')]/4(l  -  *  Vo ) 

P  =  3©To/to[l  +  3g'(l-m')/^5]/4(l-kVo) 


where  po  =  cos(^).  These  expressions  are  identical  to  the  standard  Eddington 
approximation  but  use  modified  definitions  of  the  optical  depth  and  scattering  albedo  to 
account  for  the  delta-function  component  of  the  phase  function,  defined  as 


t'  =  {l-Q}Qf)T 


(3.10a) 

(3.10b) 
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and  the  coefficients  Ci  and  C2  are  obtained  from  the  boundary  conditions  on  the  radiance, 
which  are 

=  Fq^Iq  at  the  upper  boundary 
F^  =  a^F^  at  the  ground 

The  expressions  can  be  generalized  to  a  multi-layer  situation  by  demanding  continuity  of 
the  fluxes  across  each  layer  interface  to  obtain  a  coupled  set  of  equations  for  the 
coefficients  within  each  layer  (Shettle  and  Weinman,  1970). 

A  preliminary  estimate  of  the  effect  of  dust  loading  on  the  solar  radiation  was 
obtained  from  a  single  layer  calculation  using  a  homogeneous,  monodisperse  particulate 
concentration  over  a  1km  deep  layer.  The  optical  depth  for  a  series  of  particle  sizes  and 
concentration  levels  typical  of  single  and  multiple  cloud  pedestals  is  shown  in  Table  3-1, 
and  it  is  clear  that  small  particles  can  exist  at  sufficiently  high  number  density  in  the 
pedestal  of  nuclear  clouds  to  produce  a  significant  optical  depth.  The  quantitative 
radiative  effects  are  calculated  using  the  delta-Eddington  model,  which  requires  an 
estimate  of  single  scattering  albedo  and  asymmetry  factor  for  dust.  A  typical  value  of  fiJa 
for  pedestal  dust  is  0.8,  lying  between  the  value  for  soot  particles,  which  is  around  0.4, 
and  that  for  non-carbon  dust  particles,  which  have  a  scattering  albedo  of  roughly  0.95 
(Tsay  et  al.,  1991).  The  asymmetry  factor,  g,  was  taken  to  be  0.5  for  the  dust  particles;  an 
accurate  value  is  not  available  but  the  results  are  not  very  sensitive  to  g. 

Using  £Ub  =  0.8,  g  =  0.5,  and  as  =  0,  the  values  of  transmitted  and  absorbed 
fraction  for  an  incident  vertical  solar  flux  through  a  range  of  optical  depths  covering 
Table  3-1  are  given  in  Table  3-2.  The  transmitted  fraction  is  the  fraction  of  incident  flux 
passing  through  the  dust  layer  and  reaching  the  ground,  and  the  absorbed  fraction  is  the 
amount  of  radiation  absorbed  within  the  dust  layer  which  will  cause  the  layer  to  heat  up. 
The  remaining  fraction  is  reflected  back  from  the  top  of  the  dust  layer.  Table  3-2  shows 
that  for  optical  depths  greater  than  unity,  the  transmitted  flux  is  greatly  reduced  below  the 
incident  value,  and  that  the  bulk  of  the  radiation  is  absorbed  by  the  dust  layer.  Thus,  a 
1km  layer  of  lO^im  particles  at  a  concentration  of  IQ-'^g/cc  would  virtually  extinguish  the 
radiative  flux  at  the  surface,  with  only  4%  of  the  solar  flux  being  transmitted,  while  79% 
of  the  solar  radiation  would  be  absorbed  within  the  dust  layer  itself.  These  changes  in  the 
solar  heating  effects  would  dramatically  alter  the  atmospheric  boundary  layer  structure 
and  the  vertical  mixing  rates. 
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Table  3-1.  Extinction  optical  depth  for  Ikm  dust  layer  at  uniform 

concentration. 


Concentration 

10*^  g/cc 

10’^  g/cc 

10'^  g/cc 

d  =  lOOiim 

0.0066 

0.066 

0.66 

10p.m 

0.066 

0.66 

6.6 

l|im 

0.66 

6.6 

66.6 

Table  3-2.  Transmitted  /  absorbed  fraction  for  dust  layer. 


Optical  depth 

0.066 

0.2 

0.66 

1 

2 

a 

6.6 

66.6 

Transmitted, /f 

0.98 

0.93 

0.78 

0.68 

0.43 

0.16 

0.04 

0.00 

Absorbed, /a 

0.01 

0.04 

0.14 

0.22 

0.42 

0.67 

0.79 

0.83 

3.2  DISPERSION  EFFECTS. 

In  order  to  assess  the  dynamic  effects  of  the  absorption  of  solar  radiation  in  a 
dust-laden  boimdary  layer,  we  extended  the  turbulence  model  described  by  Sykes,  Parker 
and  Henn  (1993,  Chapter  4)  for  calculating  the  evolution  of  a  spectrum  of  dust  sizes  in 
the  atmospheric  boundary  layer.  A  multi-layer  version  of  the  delta-Eddington  scheme 
was  added  to  the  one-dimensional  turbulence  closure  model,  and  the  radiative  flux 
divergence  was  included  in  the  mean  temperature  equation  as  a  local  heating  term.  The 
net  radiative  flux 

Fff= 

and  the  heating  rate  is  computed  from  the  flux  divergence  as 

_L^ 

pCp  dz 
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Figure  3-2.  Particle  size  distribution  and  particle  bin  ranges  for  the  calculation  of 
radiative  effects  on  the  DICE  sweep-up  case. 


The  model  was  initialized  with  the  pedestal  dust  distribution  from  a  high- 
resolution  DICE  calculation  of  the  sweep-up  layer  at  r=5min.  from  a  1  megaton  (MT) 
burst  at  a  scale  height  of  200  ft/kT^/^.  The  horizontally-integrated  DICE  loading  profile 
must  be  scaled  to  produce  the  desired  concentration  distribution,  and  three  arbitrary  levels 
were  chosen  to  examine  the  boundary  layer  response.  The  assumed  particle  size 
distribution  is  shown  in  Figure  3-2,  with  the  10  size  bins  indicated  by  the  vertical  Unes. 
The  dust  distribution  for  the  medium  loading  level  is  shown  in  Figure  3-3,  with  each  line 
corresponding  to  the  cumulative  mass  in  increasing  size  bins,  i.e.,  the  left-most  line 
represents  the  concentration  in  the  lowest  bin,  while  the  right-most  line  represents  the 
total  over  all  bins.  The  calculation  uses  only  the  lowest  500m  of  the  dust  profile,  as  was 
used  in  Figure  4-9  of  Sykes,  Parker  and  Henn  (1993).  The  low  and  high  loading  cases  are 
obtained  from  the  medium  loading  by  reducing  and  increasing  the  concentrations  by  a 
factor  of  10,  respectively.  The  optical  depth  is  computed  as  a  sum  over  the  various  size 
bin  contributions,  i.e., 

^  =  i_ry£eLd2  (3.11) 


20 


Figure  3-3.  Initial  concentration  profiles  from  the  pedestal  (i.e.,  z<500m)  of  the  DICE 
sweep-up  calculation.  The  medium  loading  concentrations  are  shown. 

where  a  represents  the  size  bin,  and  c  and  d  are  the  mass  concentration  and  particle 
diameter,  respectively. 

The  meteorological  conditions  can  vary  widely  as  the  time  of  release,  geostrophic 
wind,  surface  roughness,  etc.  are  changed.  We  examined  two  typical  wind  speeds, 
10ms“i  and  2ms“^  with  a  fixed  roughness  of  Im  (representing  forested  terrain)  and  an 
early  morning  and  a  midday  release,  as  described  by  Sykes,  Parker  and  Henn  (1993),. 
The  release  time  affects  the  solar  radiation,  and  we  use  a  linear  variation  of  60  between 
-90°  and  +90°  during  the  daytime  hours  between  0600  and  1800.  This  implies  a 
sinusoidal  variation  of  the  effective  downward  radiative  flux  at  the  top  of  the  model,  and 
we  take  the  solar  constant,  Fq,  to  be  lOOOWm"^  The  surface  albedo  is  arbitrarily  taken  to 
be  0.5  for  these  calculations,  implying  that  50%  of  the  incident  radiation  is  absorbed  at 
the  surface.  The  surface  albedo  can  vary  widely  according  to  the  nature  of  the  ground 
cover,  but  a  specific  value  is  required  for  the  model  calculation  and  0.5  represents  a 
moderate  albedo  for  illustrative  purposes.  As  discussed  in  section  2,  the  sensible  heat 
flux  into  the  atmosphere  is  obtained  from  a  surface  energy  balance  which  determines  the 
fraction  conducted  into  the  soil  layer  and  the  fraction  used  to  evaporate  moisture  from  the 
surface.  For  simplicity  in  the  calculations  reported  here,  we  assume  that  30%  of  the 
radiative  flux  reaching  the  surface  is  emitted  as  sensible  heat  flux  into  the  air.  This  is  a 
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crade  approximation,  but  is  sufficient  for  the  purpose  of  illustrating  the  dynamic  effects 
of  a  dust  layer  on  the  boundary  layer  turbulence. 

The  low,  medium,  and  high  dust  loadings  correspond  to  initial  optical  thicknesses  of 
about  0.6,  6,  and  60  respectively.  We  therefore  expect  the  low  dust  loading  to  show 
minimal  effect  on  the  boundary  layer,  while  the  high  loading  should  produce  a  large 
modification.  These  expectations  are  confirmed  in  the  results  for  an  0600  release  in  a 
2ms“i  geostrophic  wind  as  shown  in  Figure  3-4.  The  figure  shows  all  three  loadings,  plus 
an  imperturbed  boundary  layer  calculation  (i.e.,  no  radiative  effects  due  to  the  dust)  at 
several  times  up  to  six  hours  after  release.  The  "no  radiation"  calculation  uses  tha  same 
initial  concentration  as  the  low  loading  case,  so  that  the  small  radiative  effects  can  be 
directly  compared.  We  also  show  only  the  total  concentration  from  each  run;  the  particle 
size  distributions  are  affected  somewhat  by  the  radiative  effects,  but  the  effects  are 
similar  to  those  described  in  the  boundary  layer  studies  in  Sykes,  Parker  and  Henn 

(1993). 


Figure  3-4  shows  immediate  effects  of  dust  loading  after  1  hour.  The  low  dust 
loading  shows  only  sUght  differences  from  the  zero  absorption  case  throughout  the  six 
hour  evolution,  but  the  medium  and  high  loading  cases  are  significantly  different.  The 
low  loading  case  is  rapidly  mixed  through  the  developing  convective  boundary  layer, 
while  the  medium  and  heavy  loadings  show  a  strong  suppression  of  the  mixing  in  the 
lower  part  of  the  distribution.  The  convective  turbulence,  driven  by  the  surface  heating, 
mixes  the  pedestal  through  a  depth  of  1.5km  during  the  morning  period  for  the  zero  and 
low  loadings.  In  contrast,  the  heavier  dust  loading  prevents  the  radiation  from  reaching 
the  ground,  so  that  the  surface  heat  flux  is  much  smaller  and  the  resulting  convection  is 
much  weaker.  The  radiation  is  absorbed  in  the  upper  region  of  the  dust  layer,  especially 
for  the  high  loading,  and  lofts  a  tenuous  layer  of  material  into  the  overlying  stable  region. 
The  bulk  of  the  dust  is  apparently  relatively  undisturbed  in  the  region  near  the  surface. 

An  understanding  of  the  dynamic  effects  of  the  radiative  flux  absorption  can  be 
obtained  from  examination  of  the  turbulence  and  heating  profiles  for  the  different  cases. 
Figures  3-5  through  3-8  show  the  potential  temperature,  the  vertical  velocity  variance,  the 
sensible  heat  flux,  and  the  radiative  flux  profiles  at  the  same  four  times  as  Figure  3-4. 
The  zero  case  shows  the  typical  morning  evolution  of  the  boundary  layer;  the  shallow 
surface-based  inversion  is  eroded  quickly,  then  the  mixed  layer  grows  rapidly  through  the 
weakly  stratified  region  marking  the  extent  of  the  previous  day's  boundary  layer.  The 
surface  heating  produces  a  uniformly  mixed  potential  temperature  with  a  strong  inversion 
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Concentration  Concentration 


Figure  3-4.  Concentration  profiles  after  a  morning  release  into  a  convective  boundary 
layer  with  2ms“i  wind  speed  showing  the  radiative  effect  of  dust  loading. 
Solid  line  is  high  loading,  long  dashes  mediiun  loading,  short  dashes  low 
loading,  and  very  short  dashes  for  low  loading  case  with  no  radiative 
effects. 
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cap  and  overlying  stable  layer.  The  effect  of  the  dust  is  to  reduce  the  surface  heating,  so 
that  the  inversion  is  eroded  more  slowly,  or  not  at  all  for  the  high  loading  case.  However, 
the  dust  absorbs  the  heat  directly,  so  at  higher  altitude  the  temperatures  are  warmer  than 
the  zero  loading  case,  and  a  nearly  mixed  region  is  formed  above  the  surface  inversion. 
The  higher  loading  cases  show  a  deeper  mixed  region  due  to  the  higher  heating  through 
direct  absorption.  The  vertical  velocity  variance,  shown  in  Figure  3-6,  is  a  direct  measure 
of  the  vertical  mixing  rate  and  is  very  much  reduced  by  the  presence  of  the  dust.  The 
zero  absorption  case  gives  a  variance  of  about  Sm^s"^  in  the  middle  of  the  mixed  layer, 
while  the  high  loading  case  only  produces  a  variance  of  about  0.4m2s-2  in  the  upper  part 
of  the  mixed  region.  Although  the  dust  absorbs  more  heat,  it  is  not  absorbed  in  such  a 
way  as  to  produce  kinetic  energy  through  the  release  of  buoyant  energy. 

The  sensible  heat  flux  is  the  buoyancy  production  term  for  the  vertical  velocity 
variance,  and  the  profiles  in  Figure  3-7  show  the  standard  linear  decay  from  the  surface 
value,  with  a  small  negative  entrainment  flux  at  the  inversion,  for  the  zero  absorption  and 
low  loading  cases.  The  higher  loading  cases,  however,  show  that  the  surface  heat  flux  is 
eliminated  completely  by  the  radiative  absorption  in  the  dust  layer,  and  only  a  very  weak 
heat  flux  is  produced  by  the  absorptive  heating.  The  radiative  flux  profiles,  in  Figure  3-8, 
show  that  the  flux  is  absorbed  over  an  increasingly  deep  layer  by  the  high  and  medium 
loadings,  but  the  heating  rate  is  uniform  over  the  layer  so  that  there  is  very  little  buoyancy 
production.  The  flux  is  just  reaching  the  surface  for  the  medium  loading,  but  is 
completely  absorbed  well  above  the  ground  by  the  high  loading.  The  absorbing  layer 
warms  at  the  same  rate  throughout,  as  indicated  by  the  constant  flux  gradient,  so  there  is 
no  temperature  contrast  within  the  layer. 

It  is  interesting  to  note  the  depth  of  the  radiative  flux  absorption  in  the  medium 
and  high  loading  cases  in  Figure  3-8,  about  1800m  and  1600m  respectively  at  t=6  hours. 
The  absorption  can  be  understood  by  reference  to  the  optical  depth  profile,  shown  in 
Figure  3-9  for  the  high  loading  case  at  6  hours  after  release.  It  is  clear  that  the  linear 
absorption  profile  corresponds  to  the  optical  depth  increasing  from  zero  up  to  about  5. 
This  is  the  optical  depth  beyond  which  virtually  no  radiation  penetrates,  so  this  is  the 
layer  that  absorbs  all  the  available  energy.  This  layer  behaves  the  same  for  both  loading 
levels;  it  simply  represents  the  quantity  of  dust  needed  to  absorb  all  the  radiation.  Given 
a  large  reservoir  of  dust  in  the  pedestal,  only  the  amount  corresponding  to  an  optical 
depth  of  about  5  wiU  be  lofted  by  the  heating,  and  the  remainder  wiU  remain  in  the  stable 

inversion  near  the  surface. 
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t  =  1  hour 


t  =  2  hour 


t  =  4  hour 


t  =  6  hour 


Figure  3-5.  Potential  temperature  profiles  after  a  morning  release  into  a  convective 
boundary  layer  with  2ms“^  wind  speed  showing  the  effect  of  dust  loading. 
Units  are  °C,  and  the  line  patterns  are  as  in  Figure  3-4. 
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1  hour 


t  =  2  hour 


Velocity  variance  Velocity  variance 


t  =  4  hour  t  =  6  hour 


Velocity  variance  Velocity  variance 


Figure  3-6.  Vertical  velocity  variance  profiles  after  a  morning  release  into  a 
convective  boundary  layer  with  2ras"^  wind  speed  showing  the  effect  of 
dust  loading.  Units  are  m^s'^,  and  the  line  patterns  are  as  in  Figure  3-4. 
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t  =  1  hour 


t  =  2  hour 


Sensible  heat  flux 


Sensible  heat  flux 


t  =  4  hour  t  =  6  hour 


Figure  3-7.  Sensible  heat  flux  profiles  after  a  morning  release  into  a  convective 
boundary  layer  with  2ras~^  wind  speed  showing  the  effect  of  dust  loading. 
Units  are  °Kms-^  and  the  line  patterns  are  as  in  Figvure  3-4. 
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t  =  1  hour  i  =  2  hour 


t  =  4  hour 


t  =  6  hour 


Figure  3-8.  Net  radiative  flux  profiles  after  a  morning  release  into  a  convective 
boundary  layer  with  2ms"i  wind  speed  showing  the  effect  of  dust  loading. 
Units  are  °Kms-i ,  and  the  line  patterns  are  as  in  Figure  3-4. 
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Figure  3-9.  Profile  of  the  optical  depth  from  the  high  loading  case  at  t=6  hours. 

The  concentration  profiles  are  less  sensitive  to  radiative  effects  for  a  noon  release, 
as  shown  in  Figure  3-10.  In  this  case,  the  turbulent  eddies  are  already  in  existence  at  the 
initial  time,  and  the  removal  of  the  surface  heating  begins  a  decay  process.  However,  the 
dust  is  mixed  relatively  quickly  through  the  layer,  so  aU  four  cases  show  similar  results 
when  the  initial  scale  factor  is  taken  into  account.  The  lofting  process  is  visible  in  the 
medium  and  high  loading  cases  in  the  upper  part  of  the  profile,  and  the  mixing  has  clearly 
been  suppressed  in  the  lower  region,  but  the  differences  between  the  cases  are  much  less 
dramatic  than  the  early  morning  release. 

The  representation  of  the  detailed  effects  of  radiative  flux  modifications  due  to  the 
presence  of  dust  is  clearly  a  complex  issue.  The  delta-Eddington  model  gives  a  relatively 
simple  description  of  the  scattering  and  absorption  properties  of  a  dust  cloud,  although 
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t  =  1  hour 


L  =  2  hour 


Concentration  Concentration 


Figure  3-10.  Concentration  profiles  after  a  noon  release  into  a  convective  boundary 
layer  with  2ms“i  wind  speed  showing  the  effect  of  dust  loading. 
Concentration  units  are  IQ-’^g/cc,  and  the  line  patterns  are  as  in  Figure  3-4. 
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we  are  unable  to  provide  quantitative  comparisons  between  the  model  predictions  and 
atmospheric  observations.  The  model  results  depend  on  the  single  scattering  albedo  and 
the  asymmetry  factor,  and  these  parameters  are  not  known  accurately  for  dust  clouds. 
There  is  likely  to  be  a  range  of  values  for  different  dust  composition,  but  the  typical  value 
of  0.8  used  in  our  study  for  the  scattering  albedo  is  a  reasonable  mid-range  value.  The 
albedo  is  not  expected  to  lie  outside  the  range  between  high-carbon  values  of  0.5  and 
silicate  values  of  0.9.  Calculations  with  a  homogeneous  layer  show  that  relatively  high 
concentrations  of  small  particles  are  required  to  produce  a  strong  effect  on  the  solar 
radiative  flux.  A  concentration  of  10“^g/cc  of  10|im  particles  throughout  a  1km  deep 
layer  produces  an  optical  depth  of  6.67  and  only  4%  of  the  solar  radiation  will  penetrate 
to  the  surface.  Concentrations  of  this  magnitude  over  a  deep  layer  are  unlikely  from  a 
single  burst,  but  are  certainly  possible  for  a  multiburst  scenario.  The  radiative  flux 
model,  in  conjunction  with  a  second-order  closure  model  of  turbulent  transport  dynamics, 
has  shown  how  the  radiation  absorbed  by  the  dust  will  loft  a  fraction  of  the  cloud  The 
lofted  fraction  corresponds  to  an  optical  depth  of  about  5,  which  is  the  depth  required  to 
absorb  most  of  the  solar  flux.  No  radiation  penetrates  further  into  the  dust  cloud  and 
therefore  produces  no  heating  of  the  lowest  regions  of  the  cloud.  It  is  important  to 
remember  that  the  optical  depth  is  simply  a  geometric  quantity,  as  defined  in  this 
discussion,  and  the  actual  depth  of  penetration  will  depend  on  the  single  scattering 
albedo.  A  value  of  0.8  implies  that  80%  of  the  energy  is  scattered  by  the  particle,  so  that 
the  absorption  is  less  than  e“^’  however,  the  multiple  scattering  allows  further  absorption, 
so  the  value  is  greater  than 

The  delta-Eddington  model  could  be  implemented  in  a  full  dynamic 
meteorological  model,  but  the  implementation  in  a  Lagrangian  late-time  dispersion  model 
is  much  more  difficult.  Clearly,  the  late-time  dispersion  model  does  not  attempt  a 
detailed  description  of  dynamics,  so  we  must  consider  the  appropriate  level  of 
parameterization  for  the  radiative  effects.  There  are  two  main  phenomena  to  consider  in 
the  development  of  a  parameterization  scheme.  First,  the  effect  of  the  cloud  on  the 
surface  fluxes  and  the  resulting  boundary  layer  structure  should  be  modeled.  Second,  we 
would  like  to  represent  some  of  the  lofting  effect  of  the  solar  heating  on  the  dust  cloud. 

The  Lagrangian  framework  means  that  we  have  no  underlying  grid  on  which  to 
perform  the  vertical  flux  calculations  required  by  the  delta-Eddington  model  or  a  simpler 
estimation  scheme.  The  surface  flux  value  can  be  approximated  very  simply  using  the 
total  optical  depth  of  the  overlying  dust,  but  even  this  is  difficult  to  obtain  from  a 
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Lagrangian  collection  of  cloud  elements.  We  can  represent  the  transmitted  fraction 
crudely  as  where  the  optical  depth  is  defined  by  (3.11),  but  this  requires  a  vertical 
integration.  To  accomplish  this  efficiently,  a  horizontal  grid  must  be  defined  and  the 
Lagrangian  elements  assigned  to  each  grid  ceU.  The  grid  should  be  large  enough  to 
represent  the  boundary  layer  scale,  i.e.,  a  few  kilometers,  but  also  fine  enough  to  resolve 
the  Lagrangian  elements.  The  surface  heating  only  needs  to  be  estimated  on  the 
boundary  layer  evolution  timescale,  which  is  on  the  order  of  an  hour  or  two. 

The  dynamic  lofting  of  the  dust  is  a  more  difficult  problem  for  late-time 
dispersion  codes.  Most  codes  contain  no  dynamics,  although  SCIPUFF  (Lewellen  et  al. 
1989)  does  include  a  simplified  treatment  of  buoyant  rise.  The  description  of  radiative 
heating  is  more  complex  than  that  for  a  buoyant  source,  however,  since  the  distribution  of 
the  heating  depends  on  the  dust  distribution  itself. 
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SECTION  4 


PARTICLE  DEPOSITION 


4.1  MODEL  REVIEW. 

The  deposition  of  particles  to  the  earth’s  surface  is  important  in  determining  the 
dust  concentration  in  the  context  of  late-time  dispersion.  For  large  particles  (>30p.m), 
gravitational  settling  is  the  key  transport  mechanism  responsible  for  deposition. 
However,  for  smaller  particles,  turbulent  fluxes  and,  very  near  the  surface,  viscous 
diffusion  are  more  important  in  transporting  particles  to  the  surface.  Deposition  to  bare 
surfaces  such  as  water  or  desert  are  relatively  well  understood  and  is  mainly  a  function  of 
surface  roughness  and  wind  speed.  However,  the  situation  is  typically  much  more 
complex  since  most  of  the  earth's  land  surface  is  covered  by  some  type  of  vegetative 
canopy,  ranging  from  grasses  to  crops  to  forests.  In  this  case,  the  presence  of  a  canopy 
not  only  modifies  the  turbulent  fluxes  but  also  introduces  a  number  of  potential 
"pathways"  for  depositing  particles  to  the  canopy  elements.  Some  of  the  important 
deposition  processes  discussed  in  the  literature  are  impaction  on  surface  features  by 
turbulent  motion,  interception  by  microscale  surface  features  and  diffusion  by  Brownian 
motion.  These  wUl  be  discussed  further  in  greater  detail.  But  we  will  note  here  that  the 
relative  importance  of  these  mechanisms  depends  on  particle  size.  Brownian  diffusion, 
interception  and  impaction  are  each  of  primary  importance  for  increasingly  larger 
particles  (but  gravitational  settling  eventually  dominates  for  large  particles).  The 
Brownian  diffusion  of  small  particles  (<0.1p.m)  is  fairly  weU  understood  from  laboratory 
experiments,  but  impaction  and  interception  are  much  more  difficult  to  characterize.  See 
Slinn  (1982)  for  further  discussion.  Other  complications  such  as  phoretic  effects,  the 
presence  of  moisture  and  particle  growth  will  not  be  considered. 

Particle  deposition  is  typically  described  in  terms  of  the  deposition  velocity  , 
which  is  defined  by, 

v^  =  F,/c  (4-1) 

where  c  is  the  particle  concentration  and  Fc  is  the  total  vertical  particle  flux,  including  the 
turbulent  concentration  flux,  viscous  diffusion  and  gravitational  settling.  This  expression 
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is  derived  from  the  scalar  conservation  equation  by  assuming  stationarity  and  horizontal 
homogeneity  (Businger,  1986).  Note  that  c,  Fc  and  therefore,  vd  need  to  be  defined  at  a 
specific  reference  height,  typically  Im  to  100m  above  the  surface  or  vegetative  canopy 

height. 

Field  observations  of  deposition  velocity  for  a  variety  of  canopies  and  particle 
sizes  have  been  reviewed  by  McMahon  and  Denison  (1979),  Sehmel  (1980)  and 
Nicholson  (1988).  It  is  evident  from  these  reviews  that  there  is  a  great  deal  of 
inconsistency  between  different  experiments,  with  variations  of  two  orders  of  magnitude 
in  the  reported  values  of  Vd .  Some  of  the  inconsistencies  probably  reflect  real  variability 
in  the  environmental  conditions,  although  others  might  reflect  inadequacies  in  the 
experimental  techniques.  Businger  (1986)  presents  a  review  of  the  experimental 
techniques  used  to  measure  deposition  and  identifies  a  number  of  corrections  and 
possible  error  sources  by  examining  the  physical  processes  responsible  for  deposition. 

The  variability  in  the  observed  values  of  Vd  indicates  that  constructing  a 
theoretical  model  of  deposition  will  be  subject  to  much  uncertainty.  On  the  other  hand, 
models  which  contain  most  of  the  deposition  physics  will  be  useful  in  examining 
parametric  variations  and  uncertainties,  e.g.,  Gould  and  Davidson  (1992),  which  would 
be  difficult  to  achieve  experimentally  and  would  still  be  subject  to  uncontrolled 
environmental  variability.  We  will  not  review  the  many  theoretical  models  which  have 
been  presented  in  the  Uterature.  An  indication  of  the  various  modeling  assumptions  and 
varying  levels  of  complexity  can  be  gleaned  from  the  following  papers:  Lewellen  and 
Sheng  (1980),  Bache  (1979a,b),  Davidson,  eL  al.  (1982),  Slinn  (1982),  Wiman  and  Agren 
(1985)  and  Peters  and  Eiden  (1992). 

Deposition  velocity  is  often  modeled  in  terms  of  the  resistance  of  the  surface  layer 
to  deposition.  Resistance  is  then  defined  as  the  reciprocal  of  the  deposition  velocity, 
R  =  where  the  deposition  velocity,  v</ ,  is  redefined  as  the  excess  over  that  due  to 
gravitational  settling.  This  concept  is  convenient  since  the  total  resistance  can  be 
expressed  as  the  sum  of  the  resistances  of  various  layers  through  which  the  particles  must 
pass  on  the  way  to  the  surface  (Lewellen  and  Sheng,  1980).  For  instance,  Lewellen  and 
Sheng  consider  four  regions  defining  the  flow  near  the  surface  and  give  model 
expressions  for  their  "aerodynamic"  resistances.  The  regions  are  the  outer  boundary 
layer,  the  constant  flux  layer  (just  above  a  smooth  surface  or  vegetative  canopy),  the 
canopy  (e.g,  crops  or  forests),  and  the  viscous  sublayer  close  to  a  smooth  surface,  be  it  the 
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flat  plate  in  a  laboratory  or  the  leaf  in  a  forest  canopy.  We  consider  the  total  aerodynamic 
resistance  and  a  flnal  surface  resistance  to  actually  depositing  particles  after  they  have 
been  brought  close  to  canopy  elements  by  the  turbulent  eddies.  Thus 

vf^=R  =  Ra  +  Rd  (4.2) 

where  Ra  is  the  aerodynamic  resistance  and  is  the  final  surface  resistance.  Ra  is 
typically  set  et^ual  to  momentum  resistance  (Lewellen,  1985),  and  so  is  the  reciprocal  of 
the  "momentum  deposition  velocity": 

R-^=u»IUr  (4.3) 

where  Ur  is  the  mean  velocity  at  the  reference  height.  As  mentioned  above,  there  are 
multiple  pathways  to  surface  deposition.  The  analogy  to  resistances  in  electrical  circuits 
suggests  that  Rd  consists  of  parallel  resistances: 

Rf^  =  RB-^  +  RiN'^  +  Rm~"  (4.4) 

/?B  is  the  resistance  of  the  viscous  boundary  layer  flow  within  a  millimeter  of  the  canopy 
element  surfaces,  is  due  to  particle  interception  and  Rim  to  particle  impaction. 

Slinn  (1982)  considers  a  different,  but  nearly  equivalent,  approach  to  particle 
deposition.  He  postulates  that  deposition  velocity  is  equal  to  a  collection  efficiency 
multiplying  the  momentum  deposition  velocity, 

=  E  u*  !  Uj.  (4.5) 

With  the  same  deposition  mechanisms  considered  above,  Slinn  sets  the  total  collection 
efficiency,  E,  to  be 

£  =  1  -  (1  -  Es){\  -  Ejf,)(l  -  Ejm)  (4.6) 

where  Eg,  etc.  are  the  individual  collection  mechanism  efficiencies.  This  approach  is 
also  taken  by  Davidson,  et.  al  (1982).  They  give  an  explicit  definition  of  local  efficiency 
77  for  a  canopy  as 

_ _ particles  depositing _ 

^  particles  passing  through  area  dfh 

such  that  the  total  particle  flux  is  given  as 

F  =  ri(z)u(z)c(z)dz  (4-8) 
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where  N  is  the  number  of  canopy  elements  (assumed  cylindrical)  of  diameter  df  per  unit 
area  of  ground,  u  is  the  mean  velocity  and  h  is  the  canopy  height. 

If  the  individual  collection  efficiencies  are  defined  by 
Eg  =  Ral Rb’  ^IN  - 

it  can  be  shown  that  the  resistance  approach  and  the  collection  efficiency  approach  are 
equivalent  for  small  collection  efficiencies.  We  will  continue  the  analysis  of  the 
individual  deposition  mechanisms  in  terms  of  collection  efficiency.  It  should  be 
remembered  that  canopy  efficiencies  are,  in  principle,  height  dependent  and  that  the 
overall  collection  efficiency  will  involve  integration  through  the  depth  of  the  canopy. 


In  the  region  very  close  to  a  relatively  smooth  surface  such  as  water,  a  flat  plate  or 
canopy  elements  such  as  leaves  or  grass.  Brownian  motion  can  be  important  in 
transporting  particles  to  the  surface  (assuming  that  surface  roughness  elements  do  no 
protrude  through  the  viscous  sublayer).  This  is  particularly  so  for  very  smaU  particles 
(diameter<0.1iim)  which,  because  of  their  small  inertia,  are  able  to  follow  small  scale 
velocity  fluctuations.  Brownian  motion  results  in  the  diffusion  of  particles,  i.e.,  transport 
from  high  to  low  concentration  regions.  As  discussed  in  Lewellen  (1985),  experimental 
evidence  indicates  that  the  deposition  due  to  Brownian  motion  can  be  expressed  in  terms 
of  the  Schmidt  number,  Sc  =  v/D,  raised  to  some  powerless  than  1;  thus 

(4.10) 

Ra 

Here  v  is  the  kinematic  viscosity  of  air  and  Z)  is  the  particle  diffusion  coefficient;  values 
of  Z)  as  a  function  of  particle  radius  can  be  found  in  Friedlander  (1977).  We  follow 
Lewellen  and  Sheng  in  using  aB=0.7;  this  is  close  to  Slinn's  value  of  2/3  but  does  not 
imply  unwarranted  precision .  Then,  in  terms  of  collection  efficiency 

where  is  the  ratio  of  the  viscous  drag  to  the  total  canopy  drag  including  pressure 
forces.  This  factor  is  needed  since  the  particle  deposition  scales  only  with  the  momenmm 
transport  due  to  viscous  forces. 


Interception  is  essentially  filtration  in  that  particles  are  deposited  when  they  are 
within  a  particle  radius  of  the  collecting  elements  (the  characteristic  size  of  the  element  is 
much  greater  than  the  particle  radius).  Particle  inertia  is  ignored  so  that  the  particles  are 
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moving  with  the  fluid  arotmd  the  collecting  element.  Many  interception  models  are 
based  on  the  filtering  efficiency  of  cylindrical  fibers  in  potential  flow, 

E„=—  (4.12) 

“f 

where  a  is  the  particle  radius  and  Cf  is  the  filter  fiber  radius  (Fuchs,  1964).  However,  for 
the  more  complex  situation  in  a  canopy,  Sliim  considers  both  "small"  collectors  such  as 
vegetative  hairs  and  "large"  collectors  such  as  grass  blades  and  pine  needles  and  proposes 

'fi  ^  ^ 


^IN  - 


a 


a  +  A, 


a 


a  +  A 


LJ 


(4.13) 


where  As  is  the  characteristic  radius  of  "small"  collectors,  /  is  the  fraction  of  total 
interception  by  these  collectors  and  Ax,  is  the  characteristic  radius  of  large  collectors. 
Slinn  (1982)  shows  that  this  formulation  is  in  good  agreement  with  Chamberlain's  (1967) 
wind  tunnel  measurements  of  deposition  to  grass  for  his  chosen  parameters.  Thus,  we 
will  use  Slinn's  expression  while  recognizing  that  it  is  crude  and  requires  a  number  of 
empirical  assumptions.  Some  models  ignore  this  mechanism  completely,  e.g.,  Peters  and 
Eiden  (1992)  who  argue  that  it  is  negligible  in  their  particular  application  to  a  spruce 
forest 


Particles  are  subject  to  impaction  when  they  are  much  smaller  than  the  collection 
elements  but  large  enough  so  that  their  inertia  prevents  them  from  following  the  fluid 
flow  around  the  elements.  Collection  efficiency  by  impaction  is  a  function  of  a  Stokes 
number,  which  is  the  ratio  of  the  particle  "stop-distance"  to  a  characteristic  canopy 
element  length.  Slirm  (1982)  uses  a  definition  involving  the  friction  velocity  as  do  Peters 
and  Eiden  (1992).  This  is  in  contrast  to  Davidson  et.  al.  (1982),  who  use  the  mean  wind 
velocity.  Following  Slirm,  we  define  the  Stokes  munber 


St  = 


cAj 


cAi 


(4.14) 


where  c  is  an  empirical  constant  set  to  1  and  Tg  is  the  particle  relaxation  time.  The 
impaction  efficiency  is  given  by  Slirm  as 

.2 


^IM  - 


St^ 


i+sr 


(4.15) 


This  is  a  rather  simple  relationship  compared  to,  for  instance,  Davidson  and  Friedlander 
(1978),  who  give 
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(4.16) 


Sf 


St^  +  0.7535?^  +  2.196St  -  0.202 

based  on  a  fit  to  laboratory  filtration  data.  Another  model,  used  by  Bache  (1979b), 
Wiman  and  Agren  (1985)  and  Peters  and  Eiden  (1992),  is 


- 


St 


Yim 


(4.17) 


where  aim  and  him  are  constants  (set  to  0.8  and  2  respectively  in  Peters  and  Eiden  (1992)). 
Again,  we  follow  Slinn's  model  in  view  of  the  demonstrated  fit  with  Chamberlain's  data. 


In  situations  where  there  is  no  canopy,  e.g.,  deposition  to  water,  bare  soil  or  rock, 
some  modifications  to  the  collection  efficiency  terms  are  necessary.  Naturally,  there  is 
no  interception,  so  E/w=0.  For  Eb  and  Em,  we  use  a  generalization  of  Lewellen  and 
Sheng's  model  for  particle  deposition  to  a  smooth  flat  surface.  Thus,  we  set 


Eb  =  0.85c"®'^ 

(4.18) 

and 

r.  _  ^IM 

(4.19) 

"  1  4.  4  ' 

1  +  A/jif 

where 

A;A^  =  0.085r(l-e^-^^^') 

(4.20) 

In  this  case,  the  Stokes  number  is  given  by 

uhg 

St^ — ^ 

(4.21) 

V 


Predictions  from  this  model  compare  favorably  with  Sehmel's  (1973)  experimental  data 
on  deposition  to  a  smooth  brass  surface.  It  should  be  noted  that  in  modifying  Lewellen 
and  Sheng's  expression,  we  have  used  the  quasi-equilibrium  assumption  q  =  5.66m*  , 
where  is  twice  the  turbulent  kinetic  energy.  Also,  Lewellen  and  Sheng  suggest  that  for 
particle  diameters  larger  than  about  lOfxm,  (4.20)  should  be  divided  by  a  factor 
+  K^qXgj to  account  for  the  inability  of  the  particles  to  follow  the  turbulent 
eddies.  Here,  A  is  the  length  scale  of  the  turbulent  eddies  and^r4  =  1.  If  A  can  be 
estimated,  we  suggest  including  this  factor,  but  it  is  usually  very  close  to  unity  and  can 
often  be  ignored. 
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A  further  complicating  phenomenon  is  particle  rebound,  which  reduces  the 
deposition  efficiency.  It  is  quite  dependent  on  surface  characteristics  and  is  not  well 
understood  (Slinn  1982).  Slinn  defines  a  rebound  fraction  which  multiplies  the  total 
deposition  efficiency.  He  suggests  a  simple  function  given  by 

=  (4.22) 


where  bis  a.  constant  Since  the  understanding  of  rebound  is  quite  poor,  we  typically  set 
b=0,  except  in  the  case  of  Chamberlain's  (1967)  wind  tuimel  data,  described  below,  where 
b=2  as  in  Sliim  (1982). 


4.2  MODEL  DEVELOPMENT. 

To  compute  the  actual  deposition  flux,  the  collection  models  just  discussed  must 
be  coupled  with  models  or  measurements  of  momentum  flux  within  the  canopy.  Three 
approaches  commonly  taken  are 

1)  assiune  shapes  for  the  velocity  and  eddy  diffusivity  profiles,  e.g., 

Slinn  (1982)  and  Davidson,  et  al.  (1982), 

2)  use  measurements  of  velocity  and  make  assumptions  about  the 
diffusivity,  e.g.,  Bache  (1984)  and  Peters  and  Eiden  (1992)  and 

3)  compute  the  velocity  and  eddy  diffusivity  given  canopy 
characteristics,  e.g.,  Lewellen  and  Sheng  (1980). 

The  first  approach  is  a  kind  of  bulk  parameterization  of  the  flow  through  the  canopy 
while  the  second  is  obviously  only  feasible  if  experimental  data  is  available.  For  the 
typical  situation  where  the  complete  deposition  model  must  include  some  prediction  of 
the  flow  through  the  canopy,  we  consider  three  types  of  models  of  increasing  complexity: 
a  simple  bulk  parameterization  of  the  canopy,  a  one  dimensional  model  (horizontally 
averaged)  of  the  canopy  which  resolves  the  vertical  profiles  of  mean  velocity  and 
turbulence,  and  a  three-dimensional  model,  which  is  basically  a  number  of  one¬ 
dimensional  models  driven  by  the  three-dimensional  velocity  field  above  the  canopy. 
They  all  include  in  some  way  the  deposition  mechanisms  described  above. 
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4.2.1  Bulk  Canopy  Model. 


The  bulk  canopy  model  is  the  simplest  and  gives  deposition  velocity  as  a  function 
of  total  surface  momentum  flux.  As  mentioned  above,  this  type  of  model  involves 
assumptions  about  the  velocity  profile  and  turbulent  diffusivity  in  the  canopy.  For 
example,  Slum's  (1982)  model  requires  the  specification  of  parameters  such  as  shape 
factors,  the  ratio  of  canopy-top  velocity  to  the  velocity  at  the  reference  height  and  zero- 
plane  displacement  However,  given  the  capability  of  one-dimensional  second-order 
turbulence  closure  models,  such  as  that  to  be  described  shortly,  to  compute  these  profiles 
without  making  such  gross  assumptions,  we  suggest  that  models  such  as  Slum's  be  further 
simplified,  making  them  appropriate  for  use  in  large-area  (e.g.  continental  scale)  three- 
dimensional  models.  For  instance,  if  we  make  the  assumption  that  E  is  small  (say  0.1  or 
less),  then  Slum's  model  can  be  simplified  to  the  form 


Vd 


E 

{l-P)E+P 


ul 


— +gT, 


M. 


(4.23) 


where  )3  is  a  parameter  which  in  general  will  vary  for  different  canopies  and  velocity 
profiles.  This  is  effectively  defining  an  overall  collection  efficiency  for  the  canopy. 
Slinn  gives  some  data  which  indicate  that  the  range  of  p  is  approximately  0.1  to  0.3.  Fits 
to  the  data  of  Chamberlain  (1967)  using  (4.23)  with  ^=0.16,  along  with  the  full  Slinn 
model  (see  SUnn  (1982)  for  the  full  expression),  are  shown  in  Figure  4-1.  It  is  seen  that 
the  agreement  between  the  simple  model  and  the  full  model  is  quite  good;  both  models 
are  in  good  agreement  with  the  data  as  well.  The  results  are  not  very  sensitive  to  changes 
of  ±0.05  in  p.  It  remains  to  be  determined  if  this  simple  model  can  be  applied  to  other 
canopies  and  under  a  variety  of  atmospheric  conditions. 


4.2.2  Canopy  Profile  Model. 

The  horizontally-averaged  canopy  model  solves  for  the  vertical  profiles  of  mean 
velocity,  particle  concentration  (or  some  other  scalar),  temperature,  humidity  and 
turbulent  fluctuations  using  the  full  second-order  turbulence  closure  scheme  of  Lewellen 
(1977).  The  modifications  required  to  model  the  flow  in  a  canopy  are  discussed  in 
Lewellen  and  Sheng  (1980)  and  will  only  be  briefly  summarized  here.  The  total  drag 
force  is  modeled  as 

Di  =  1 
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particle  timescale  (sec),  v^/g 

Figure  4-1.  Comparison  between  Chamberlain  (1967)  wind  tunnel  particle  deposition 
data  (symbols),  Slinn  (1982)  model  (solid  line),  and  bulk  canopy  model 
(4.23)  (dashed  Une). 


where  Cf  and  Cp  are  the  skin  friction  and  pressure  drag  coefficients,  respectively.  Ay/  is 
the  total  wetted  area  (per  volume),  is  the  frontal  area  (per  volume),  is  twice  the 
turbulent  kinetic  energy  and  w,-  denotes  the  mean  velocity  vector.  Cj  is  modeled  as 


Cf=Ci 


UaJ 


(4.25) 


where  A  is  the  turbulence  dissipation  length  scale  and  c\  is  a  constant  set  to  1. 


There  are  sink  terms  in  the  turbulent  energy,  temperature,  humidity  and  particle 
concentration  equations  which  are  modeled  simply  by  scaling  the  friction  drag  term  by  an 
appropriate  function  of  the  Prandtl  number.  The  Reynolds  stress  equations  are  also 
modified  by  adding  a  source  term  to  account  for  the  creation  of  wake  turbulence  due  to 
pressure  drag  and  a  .sink  term  which  models  the  dissipation  of  turbulent  fluctuations  by 
skin  friction.  Further  details  can  be  found  in  Lewellen  and  Sheng  (1980). 


The  one-dimensional  particle  deposition  model  as  given  here  uses  much  of  Slinn  s 
formulations  in  place  of  those  given  in  Lewellen  and  Sheng,  although  it  should  be 
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emphasized  that  the  full  particle  concentration  profile  through  the  canopy  is  being  solved 
for.  Thus,  we  are  merely  modifying  the  sink  term  in  the  conservation  equation  of 
Lewellen  and  Sheng,  which  is  written  as 

Here,  cd  is  the  total  drag  coefficient  such  that  the  momenmm  flux  due  to  the  canopy  is 
—Cdu,  E  is  the  total  collection  efficiency  as  discussed  above  and  K  is  the  eddy  diffusivity. 
It  should  be  noted  that  collection  efficiencies  are  local,  i.e.,  height  dependent,  as  in  the 
model  of  Davidson,  et.  al.  (1982)  and  Peters  and  Eiden  (1992).  In  this  context  we  define 
the  local  Stokes  number  for  use  in  to  be 
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Davidson,  et.  al.  (1982)  use  the  local  mean  velocity  while  Peters  and  Eiden  (1992)  define 
a  "local"  M*. 

The  results  of  the  one-dimensional  model  for  the  grass  canopy  of  Chamberlain 
(1967)  are  shown  in  Figure  4-2,  along  with  Chamberlain’s  data.  The  vertical  extent  of  the 
computational  domain  is  twice  the  canopy  height  of  6cm  and  is  resolved  by  18  grid 
points.  The  assumed  distribution  of  is  shown  in  Figure  4-3.  The  ratio  A^/A^is 
assumed  to  be  3.  The  agreement  with  the  data  is  good,  indicating  that  the  model  is 
correctly  including  much  of  the  physical  processes  occurring  in  the  canopy.  The  simpler 
bulk  canopy  model  matches  the  data  somewhat  better,  but  it  should  be  noted  that  the 
adjustable  parameters  in  Slinn's  model  and  the  related  parameter  ^  were  tuned  to  agree 
with  the  data.  The  extension  to  other  canopies  and/or  different  wind  profiles  must  be 
made  with  caution.  The  one-dimensional  model,  on  the  other  hand,  is  more  general  in 
principle  since  it  has  fewer  adjustable  parameters  (or  at  least  fewer  "critical  parameters). 
The  second-order  closure  turbulence  model  parameters  are  determined  from  fits  to 
standard  laboratory  flows.  The  canopy  terms  require  an  estimate  of  the  pressure  drag 
coefficient  (and  canopy  area  density  which,  ideally,  is  measured  in-situ).  Although  the 
drag  coefficient  is  undoubtedly  not  the  same  for  all  canopies,  Lewellen  and  Sheng  (1980) 
show  that  Cp=0.16  gives  reasonable  results  for  both  com  and  forest  canopies.  (The 
foregoing  applies  principally  to  modeling  mean  velocity.  The  deposition  mechanisms  are 
much  less  well  understood  and  are  subject  to  greater  uncertainty.)  Thus  it  can  be  applied 
to  other  canopies  with  some  degree  of  confidence,  as  we  now  show. 
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particle  timescale  (sec), 

Figure  4-2.  Comparison  between  Chamberlain  (1967)  wind  tuimel  particle  deposition 
data  (symbols)  and  explicit  one-dimensional  canopy  model  (4.26)  (solid 
line). 


Figure  4-3.  Assumed  grass  leaf  area  density  profile  for  comparison  with  Chamberlain 
laboratory  data,  made  non-dimensional  by  canopy  height,  h. 
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Figure  4-4  Assumed  leaf  area  density  profile  for  three  forest  canopy  types,  made  non- 
dimensional  by  canopy  height,  k 


Although  we  are  not  aware  of  any  other  deposition  measurements  comparable  to 
Chamberlain's  there  are  quite  a  few  measurements  of  mean  velocity  and  turbulence 
statistics  presented  in  the  literature.  We  have  chosen  to  compare  with  Amiro's  (1990)  data 
for  three  boreal  forest  canopies,  pine,  spruce  and  aspen.  Based  on  Amiro’s  measurements 
of  canopy  foliage,  we  constructed  idealized  profiles  of  the  canopy  area  densities  which 
are  shown  in  Figure  4-4.  In  aU  cases  it  is  assumed  that  A^/Ay  =  4.  The  computational 
domain  is  the  same  as  for  the  grass  calculation  (suitably  scaled  by  canopy  height).  It 
should  be  noted  that  the  calculations  were  made  with  a  different  turbulence  dissipation 
length  scale  equation  than  that  of  Lewellen  and  Sheng  (1980).  An  analysis  of  the  data 
indicated  that  the  length  scale  should  be  modeled  as 

^  ^  0.05(LA/-H£)  (4  28) 

CpAfiz) 


where  eis  a  small  number  set  to  0.1,  Cp  is  the  pressure  drag  coefficient,  and  LAI,  the  leaf 
area  index,  is  defined  as 


(4.29) 


This  formulation  differs  from  that  given  in  Lewellen  and  Sheng,  where 
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a 


(4.30) 


and  0=0.1.  In  both  cases,  A  is  constrained  by  the  condition 


\dz 


<0.65 


as  suggested  by  WUson  and  Shaw  (1977).  The  original  length  scale  (4.30)  was  used  in 
Lewellen  and  Sheng's  calculation  for  a  com  canopy  and  found  to  be  quite  satisfactory. 
However,  the  formulation  (4.28),  based  on  Amiro's  forest  data,  is  not  entirely  satisfactory 
since  it  relates  the  local  scale  to  a  gross  canopy  measure,  LAI.  It  is  equivalent  to  setting 
the  constant  a  in  the  original  formulation  to  0.05*LAI.  (The  values  of  LAI  for  the  pine, 
spruce  and  aspen  canopies  are  2,  4  and  10  ,  respectively  (Amiro,  1990).)  It  was  found, 
though,  that  the  com  canopy  model  calculations  with  (4.28)  did  not  differ  much  from 
those  with  the  original  expression.  In  contrast,  calculations  of  the  forest  canopy  with 
(4.30)  were  quite  different  and  showed  poor  agreement  with  the  data. 


Figure  4-5  shows  a  comparison  between  the  one-dimensional  model  and  Amiro's 
measurements  of  mean  velocity  profiles.  Despite  the  fact  that  there  was  some  tuning  in 
the  length  scale  equation  to  match  the  experiment,  it  is  encouraging  that  the  model 
predicts  the  mean  velocity  so  well  for  different  types  of  forest  canopies.  The  vertical 
shear  stress  profiles  in  Figure  4-6  are  also  in  good  agreement  with  the  data.  In  particular, 
the  model  predicts  nearly  uniform  velocity  and  small  shear  stress  seen  in  the  bottom  half 
of  the  three  canopies  and  even  suggests  the  presence  of  a  local  velocity  maximum  near 
the  surface  for  the  spmce  and  pine  forests.  Regions  of  small  shear  or  local  velocity 
maxima  have  been  frequently  observed  in  canopies  with  reduced  foliage  density  near  the 
surface  (Shaw,  1977).  Thus,  although  this  effect  will  depend  on  our  assumed  canopy  area 
distribution,  it  is  nonetheless  a  critical  result  since  it  cannot  be  predicted  by  simple 
mixing-length  models  (Shaw,  1977)  and  is  certainly  difficult  to  fit  into  a  framework 
requiring  velocity  profile  assumptions  as  in  Slirm  (1982). 

Figure  4-7  shows  model  predictions  and  experimental  measurements  of  the 
normal  velocity  turbulent  flucmations  in  the  three  forest  canopies  calculation.  The 
vertical  velocity  variance  shows  good  agreement  with  the  data.  The  horizontal 
components  show  somewhat  larger  discrepancies,  but  are  still  in  reasonable  agreement. 
The  calculation  of  higher-order  turbulence  statistics  is  generally  more  sensitive  to 
modeling  assumptions  than  the  mean  velocity.  It  should  be  noted,  however,  that  it  is  also 
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Figure  4-5  Profiles  of  mean  velocity  for  three  forest  types  normaHzed  by  the  velocity 
above  the  canopy  from  Amiro  (1990)  compared  with  one-dimensional 

canopy  model  and  LES. 


Figure  4-6.  Profiles  of  turbulent  shear  stress  for  three  forest  types  normalized  by  the 
value  above  the  canopy  from  Amiro  (1990)  compared  with  one¬ 
dimensional  canopy  model  and  LES. 
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Figure  4-7.  Profiles  of  turbulent  velocity  fluctuations  for  three  forest  types  normalized 
by  the  value  above  the  canopy  from  Amiro  (1990)  compared  with  one¬ 
dimensional  canopy  model  and  LES. 
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difficult  to  obtain  reliable  measurements  as  weU  and  this  is  reflected  in  the  large 
uncertainties  associated  with  the  data;  typical  standard  deviations  are  about  0.1  to  0.2.  It 
is  possible  that  canopy  inhomogeneities,  non-stationarity  and/or  stability  effects  may 
result  in  measurement  uncertainty. 

4.2.3  LES  Canopy  Model. 

The  canopy  equations  have  also  been  implemented  in  a  three-dimensional  model, 
where  the  presence  of  the  canopy  effectively  modifies  the  surface  boundary  conditions. 
The  results  from  large-eddy  simulations  for  the  pine  forest  are  also  shown  in  Figures  4-5, 
4-6  and  4-7.  Here,  the  canopy  is  resolved  by  about  10  points  in  the  vertical.  The 
horizontal  grid  spacing  is  about  eight  times  the  canopy  height.  It  can  be  seen  that  the 
one-dimensional  and  three-dimensional  calculations  are  fairly  close.  There  are  some 
differences  in  Cy  and  Oiv  in  the  lower  half  of  the  canopy,  where  the  LES  predicts  greater 
transverse  fluctuation  and  reduced  vertical  fluctuations.  It  should  be  noted  that  we  are 
not  resolving  small  scale  eddies,  e.g.,  on  the  scale  of  the  canopy  height  or  less,  but  the 
large  scale  eddies  from  the  atmosphere  above  the  canopy  effectively  impose  local  mean 
wind  boundary  conditions  at  the  canopy  top. 


4.3  DEPOSITION  UNDER  CONVECTIVE  CONDITIONS. 

LES  with  an  explicit,  vertically  resolved  canopy  is  a  computationally  expensive 
model.  Late-time  atmospheric  dispersion  calculations  will  be  performed  with  a  bulk 
canopy  model  as  described  above,  and  only  mean  velocities  will  be  available  (as  opposed 
to  the  explicit  turbulence  of  LES).  The  bulk  deposition  model  requires,  at  a  minimum, 
the  surface  momentum  flux,  ,  which  will  be  determined  by  the  mean  wind  above  the 
canopy  and  the  surface  roughness  height  as  a  function  of  canopy  height  and  foliage 
density,  etc.  However,  under  free  convection  conditions,  vamshes  since  the  mean 
wind  is  zero  and  so  no  deposition  will  be  predicted.  This  is  obviously  not  the  case  in 
reality  since  individual  convective  eddies  will  transport  momentum  (and  particles)  to  the 
surface.  It  is  possible,  however,  to  define  the  average  local  friction  speed  ,  v*,  as  a 
measure  of  the  surface  momentum  flux.  As  part  of  this  study,  we  have  used  the  LES 
technique  to  calculate  the  dependence  of  v.  on  the  surface  roughness  and  convective 
velocity  scale  w*  =  where  H  is  the  surface  temperature  flux,  Tq  is  the 

reference  temperature  (300°K),  and  zi  is  the  convective  boundary  layer  height.  This  work 
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Figure  4-8.  Variation  of  the  magnitude  of  the  average  surface  friction  velocity  with 
roughness  length  for  free  convection  conditions. 

has  been  published  in  the  peer-reviewed  literature  (Sykes,  Henn  and  Lewellen,  1993).  By 
postulating  a  balance  between  inertial  terms  due  to  the  large  scale  eddies,  which  scale 
with  zi  and  the  vertical  stress  gradients  across  a  thin  surface  layer  (typically  on  the  order 
of  10“^  Zj),  a  simple  algebraic  scaling  relationship  is  obtained.  The  required  constants  are 
determined  using  LES  with  explicit  ensemble  averages  of  instantaneous  friction  velocity 
being  calculated.  The  resulting  relationship  is 

21n-^^ — lna2~  =aik  (4.31) 

W*  W*  Zi  j 

where  Zq  is  the  roughness  height  and  k  is  von  Karman's  constant  The  empirical 
constants,  =  1  and  a2  =  2,  are  determined  from  a  least-squares  fit  to  the  LES  ensemble 
averages.  The  resulting  curve  is  plotted  in  Figure  4-8.  This  curve  can  be  used  to 
determine  a  mean  friction  velocity  under  free  convection  conditions,  which  can  then  be 
used  in  a  simple  deposition  model.  A  simple  curve  fit  to  this  transcendental  equation  is 

given  by 
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(4.32) 


Probability  distributions  of  the  instantaneous  friction  velocity  normalized  by  the 
ensemble  mean  value  are  shown  in  Figure  4-9(a).  The  shapes  are  very  similar  for 
different  roughness  lengths  and  are  nearly  symmetric  about  their  means.  It  appears  that 
distributions  are  very  nearly  Gaussian,  with  the  standard  deviations  being  approximately 
one  third  of  the  mean  (but  decreasing  slowly  with  increasing  roughness).  Figure  4-9(b) 
shows  probability  distributions  for  cases  with  mean  geostrophic  wind,  but  which  are  still 
dominated  by  convective  eddies.  It  is  evident  that  the  distributions  become  narrower 
about  the  mean  as  the  geostrophic  wind  increases.  It  is  also  interesting  to  note  that  these 
simulations  indicate  that  for  even  small  mean  winds,  it*  is  very  close  to  v*  Thus  it  is 
most  likely  that  the  free  convection  v*  can  be  used  as  a  minimum  on  the  friction 

velocity: 

It*  =  max(v*,u^“^)  (4.33) 

where  'P  is  the  normalized  velocity  profile  given  by 

'P  =  r^ln— ("^-34) 
.  ^ 

and  is  a  stability  correction  to  the  neutral  logarithmic  profile  based  on  Monin- 
Obukov  similarity  theory  (Businger,  1973). 


%• 
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Figure  4-9.  Probability  density  function  for  the  instantaneous  friction  velocity. 

(a)  variation  with  roughness  length  for  free  convection  conditions; 

(b)  variation  with  geostrophic  wind  speed,  dimensionless  roughness  of 
10^  except  where  indicated. 


51 


SECTION  5 

TERRAIN  EFFECTS 


5.1  GENERAL  REMARKS. 

The  PBL  representation  discussed  in  Section  2  is  appropriate  for  a  homogeneous, 
flat  surface.  The  dynamics  of  the  planetary  boundary  layer  under  these  conditions  are 
generally  well  understood,  at  least  for  neutral  and  unstable  conditions,  and  the  mean  flow, 
turbulence,  and  dispersion  characteristics  can  be  determined  in  terms  of  a  few  important 
PBL  parameters.  The  turbulence  generation  mechanisms,  namely  shear  and  buoyancy, 
are  adequately  represented  by  current  parameterization  methods  using  boundary  layer 
quantities  such  as  the  surface  friction  velocity  and  the  convective  velocity  scale, 
(Deardorff,  1972).  Research  continues  to  reveal  details  of  the  turbulence  structure,  but 
the  main  features  of  the  flow  are  now  well  accepted. 

The  actual  situation  in  the  atmosphere  is  rarely  ideal,  however.  The  surface  is  not 
usually  homogeneous,  but  contains  variations  in  character  and  elevation  on  all  length 
scales.  The  distinction  between  the  surface  roughness  and  the  surface  elevation 
variations  is  somewhat  arbitrary,  but  the  roughness  is  usually  taken  to  refer  to  complex 
surface  features  much  smaller  than  the  boundary  layer  depth.  The  local  surface 
roughness  is  therefore  determined  by  the  surface  coverage,  e.g.,  type  of  vegetation, 
density  of  buildings,  water  surface,  etc.  The  surface  coverage  is  generally 
inhomogeneous  on  scales  of  interest  for  meteorological  prediction  models,  both  synoptic 
and  mesoscale,  so  representations  are  required  for  larger  scale  area-averages.  Recent 
efforts  to  describe  the  effects  of  non-uniformity  in  local  surface  roughness  have  been 
made  by  Mason  (1988),  Vihma  and  Savijarvi  (1991),  and  Wood  and  Mason  (1991). 


The  effects  of  terrain  elevation  variations  have  often  been  represented  using  the 
same  framework  of  surface  roughness  specification  as  the  smaller  scale  features.  Many 
studies  of  neutral  boundary  layer  flow  over  terrain  features  have  been  conducted,  and 
Taylor  et  al.  (1989)  present  a  synthesis  of  several  numerical  modeling  approaches  to 
determine  a  parameterization  of  the  drag  forces.  AU  of  the  models  employ  empirical 
closure  techniques  for  the  turbulent  flow,  although  there  is  some  insensitivity  to  the  level 
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of  closure  in  the  drag  results.  There  are  also  a  number  of  theoretical  analyses  of  neutral 
flow  over  idealized  terrain,  e.g.,  Jackson  and  Hunt  (1975),  Sykes  (1980),  Hunt  et  al 
(1988),  which  have  helped  to  delineate  the  important  phenomena  in  the  flow.  In  general, 
however,  there  are  no  reliable  parameterizations  of  the  boundary  layer  over  complex 
terrain.  We  attempt  to  improve  the  situation  by  means  of  numerical  simulation  of  the 
turbulent  processes  over  hilly  terrain. 

Our  approach  to  the  problem  is  to  use  the  best  available  turbulence  modeling 
technique  to  calculate  boundary  layer  flow  over  idealized  terrain  shapes  for  a  range  of 
meteorological  conditions.  The  numerical  technique  is  Large  Eddy  Simulation  (LES), 
which  involves  a  three-dimensional,  time-dependent  integration  of  the  equations  of 
motion,  using  a  finite-difference  grid  of  sufficient  resolution  that  the  energy-containmg 
turbulent  eddies  can  be  represented  explicitly  (Rogallo  and  Moin,  1984).  This  allows  a 
direct  computation  of  the  dynamics  of  the  large  eddies  and  avoids  any  empirical  closure 
assumption.  It  is  impossible  to  resolve  the  entire  spectrum  of  motions  m  the  atmosphere, 
so  the  T  FS  requires  a  subgrid  closure  to  represent  the  effects  of  eddies  too  small  to  appear 
on  the  finite-difference  grid.  However,  the  subgrid  closure  is  only  important  for  the 
small-scale  features  and  becomes  less  important  as  the  grid  is  refined. 

We  use  LES  to  obtain  detailed  results  for  simple  terrain  and  develop 
parameterizations  based  on  these  idealized  cases.  We  shall  examine  the  effects  of  terrain 
slope,  wavelength,  height,  surface  roughness,  and  atmospheric  stability.  The 
parameterizations  will  be  written  in  terms  of  general  terrain  characteristics,  so  that  the 
results  can  be  extended  to  more  complex  topography.  We  shall  then  test  the 
parameterizations  by  actual  computation  with  different  terrain  shapes. 


The  technique  of  Large-Eddy  Simulation  allows  a  more  reUable  representation  of 
the  turbulent  processes  than  one-point  closure  models,  provided  that  the  dominant  eddies 
can  be  resolved  by  the  calculation.  This  requirement  is  much  easier  to  achieve  in  the  case 
of  convective  boundary  layer  flows,  where  the  turbulence  generation  mechanism 
produces  large-scale  eddies  directly.  Studies  of  free  convection  flows  (Mason,  1989; 
Schmidt  and  Schumann,  1989;  Sykes  and  Henn,  1989)  have  demonstrated  the  capability 
of  LES  for  these  flows.  The  neutral  flow  situation  is  more  difficult  to  calculate  with  LES, 
since  the  shear  generation  mechanisms  are  concentrated  in  the  small  scales  near  the 
surface  where  the  shear  is  highest.  LES  has  difficulty  representing  this  region  reUably 
(Mason  and  Thompson,  1992).  In  contrast  to  the  LES  limitations,  the  theoretical 
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approaches  used  for  neutral  boundary  layer  flow  are  not  easily  extended  to  the  convective 
case.  The  theories  rely  on  a  linearization  about  the  undisturbed  boundary  layer  flow,  but 
the  effects  of  surface  heating  over  a  sloping  surface  can  produce  vertical  motions  in 
convective  cells  with  preferential  aUgnment  The  convective  flow  situation  thus  presents 
a  suitable  candidate  for  LES  investigation. 

Free  convection  over  a  wavy  surface  has  been  studied  recently  by  Krettenauer  and 
Schumann  (1992)  and  Walko  et  al.  (1992)  using  LES.  The  conclusions  from  these 
studies  are  generally  concerned  with  the  modification  of  the  convective  cell  structure  by 
the  terrain  slopes,  which  are  two-dimensional  with  no  variation  in  the  y-direction.  In  the 
absence  of  any  mean  wind,  the  terrain  can  preferentially  lock  the  updrafts  in  the  form  of 
rolls  aligned  with  the  surface  terrain  in  the  lower  part  of  the  mixed  layer,  although  the 
general  turbulence  profiles  do  not  appear  markedly  different  from  convection  over  a 
uniform  surface. 

In  the  present  study,  we  are  interested  in  the  effects  of  two-dimensional  terrain 
variations  on  convective  boundary  layer  flow  with  a  mean  geostrophic  wind.  We  shall 
restrict  attention  to  a  simple  sinusoidal  terrain  shape  and  attempt  to  determine  the 
variations  with  terrain  wavelength  and  slope.  The  principal  features  of  mterest  are  the 
mean  flow  and  turbulence  modifications.  The  integrated  momentum  deficit  (relative  to 
the  geostrophic  wind)  must  balance  the  surface  forces  in  the  steady  state,  and  this  balance 
will  be  used  to  check  the  steadiness  of  the  solutions.  The  dependence  of  the  drag  forces 
on  terrain  and  surface  parameters  will  be  determined  from  the  LES  results. 


5.2  NUMERICAL  MODEL  DESCRIPTION. 


The  basic  model  used  for  this  study  is  an  extension  of  the  Cartesian  model  of 
Sykes  and  Henn  (1989).  The  model  is  a  second-order  accurate,  finite  difference 
representation  of  the  filtered  equations  of  motion  for  an  incompressible,  Boussinesq  flmd. 
A  terrain-following  coordinate  transformation  was  introduced  using  the  techniques  of 
Clark  (1977).  The  vertical  coordinate  used  in  the  model,  C .  is  defined  as 

r  =  (5.1) 

^  J 

where 
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(5.2) 


Here  h{x,y)  is  the  local  terrain  elevation,  z  is  the  Cartesian  vertical  coordinate,  and  D  is 
the  depth  of  the  model  domain.  Following  Clark,  we  retain  the  Cartesian  velocity 
components  for  the  momentum  and  Reynolds  stress  representations  and  define  a 
transformed  vertical  component,  (O,  which  simplifies  the  advective  transport 
computation.  Thus,  (u.v.w)  are  the  velocity  components  in  the  (x,y,z)  directions 
respectively,  and  (O  is  defined  as 

CD  =  w-  GjcU  -  GyV 


where 


(5.4) 


The  momentum  equations  can  be  written 

—  =  -—{Jp)-fj(v-Vg)  +  -^{Jtn)+-^{J'^i2)+-^{'^i3+Gx'^n  +  Gyti2-GxP) 

Dt  dx  ^  dx  dy  d^ 

(5.5a) 

^  =  -—{Jp)  +  fjiu  -  M J  +  4:('^^12)  +  T-(-^'^22)  +  +  <^y^22  "  ^yP) 

Dt  dy  dx  dy 

(5.5b) 

(5.5c) 


where  the  material  derivative  is  defined  as 

:££  =  |.(7(P)+^(7m(P)+|-(7v9))+^(©(P)  (5-6) 

Dt  dt  dx  dy  dig 

In  the  above  equations,  g  is  the  gravitational  acceleration,  0  is  the  potential  temperature. 
To  is  the  Boussinesq  reference  temperature,  p  is  the  dynamic  pressure,  and  Ty-  is  the 
subgrid  stress  tensor.  The  geostrophic  wind  is  (m^,v^)  and  /  is  the  Coriolis  parameter, 
which  we  take  to  be  IQ-V^  in  all  the  calculations  reported  below.  The  mass 
conservation  equation  is  simply 


^  tT  \  ^  /  7  \ .  n 

-(7.)+-(yv)^-  =  0 


(5.7) 
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The  potential  temperature  equation  is 

Dt  <2)c  dy  oQ 


(5.8) 


where  Hi  is  the  subgrid  heat  flux. 


The  subgrid  turbulence  model  uses  a  turbulent  kinetic  energy  transport  equation 
based  on  the  second-order  closure  model  of  Lewellen  (1977).  The  subgrid  velocity 
variance,  ,  is  obtained  from  the  following  conservation  equation 


££. 

Dt 


axj 


dx 


dy 


(5.9) 


where  F;  is  the  subgrid  flux,  and  A  is  the  subgrid  turbulence  length  scale.  The  first 
two  terms  on  the  right  hand  side  represent  shear  and  buoyancy  production  respectively, 
and  the  last  term  is  the  dissipation  rate.  The  empirical  turbulence  model  constant,  b  ,  is 
0.125  (Lewellen ,  1977). 


Subgrid  fluxes  are  modeled  as 


similarly  for  F,-,  and 


Hi=K 


it 

dxi 


(5.10) 


%=v 


^Xj  dXi 


(5.11) 


V  V  '  y 

for  the  subgrid  stresses.  The  subgrid  diffusivities  are  obtained  from  the  equilibrium 
second-order  closure  results, 

K  =  SjjqA,  V  =  S„qA  (5.12) 

with  the  stability-dependent  coefficients,  Sh  and  Sm,  are  given  by 

(5.13a) 


Sfj  — 


^  A  +  Ri{2  +  ybs) 
(A^+{A/bs-l)Ri) 
A  +  Ri 


5m  - 


'H 


(5.13b) 
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(5.14) 


The  subgrid  turbulent  Richardson  number  is  defined  as 


Ri  = 


g  dd 


Finally,  the  subgrid  length  scale  is  defined  as 


An^  A 


2 

-max 


(5.15) 


where  Aq  =  mm{ocn,ql2N),  n  is  the  normal  distance  from  the  lower  boundary,  z  =  Kx,y). 
N  is  the  Brunt- Vaisala  frequency  and  the  turbulence  constants  a=  0.65,  s=  1.8  and  A=0.75 
are  chosen  to  match  the  log-layer  profile  with  a  von  Karman's  constant  of  0.4.  This 
formulation  follows  the  philosophy  of  Mason  and  Callen  (1986)  in  distinguishing 
between  the  filter  scale,  A  ,  and  the  numerical  grid  scale.  The  specification  of  Amax  is 
problem-dependent  and  will  be  given  below. 


The  finite-difference  representations  are  similar  to  those  used  in  the  Cartesian 
model  of  Sykes  and  Henn  (1989).  Temporal  differencing  uses  the  second-order  leapfrog 
scheme  with  a  small  amount  of  smoothing,  usually  1%,  to  couple  the  two  time  levels  and 
prevent  time-splitting.  The  advective  terms  use  the  'absolutely-conserving'  scheme  of 
Piacsek  and  Williams  (1970)  to  maintain  conservation  of  momentum  and  energy. 
Diffusion  terms  use  central  differencing  with  the  du  Fort-Frankel  approximation  to 
provide  a  stable  integration. 


The  elliptic  equation  for  the  pressure  field  is  obtained  by  application  of  the  fimte- 
difference  divergence  operator  to  the  partially  advanced  velocity  field,  i.e.,  including  all 
terms  except  the  pressure  gradient.  The  pressure  field  is  then  computed  from  the  non- 
separable  equation  to  maintain  zero  divergence  of  the  velocity  field.  The  presence  of  the 
terrain  terms  prevents  the  use  of  the  orthogonal  decomposition  methods  of  the  Cartesian 
model,  so  an  iterative  scheme  is  employed.  The  pressure  equation  is  written  in  the  form 

Lo{p)  =  N{p)  +  S  (5-16) 

where  Lq  represents  the  Cartesian  finite-difference  operator  (i.e.,  the  case  h  =0),  N 
represents  the  terrain  transform  terms,  and  S  is  the  divergence  source  term.  The  solution 
is  obtained  using  the  Cartesian  solver  for  each  iteration  of  the  system 

+  S  (5-17) 
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using  the  previous  timestep  value  as  the  first  guess  for  The  iteration  is  continued 
until  the  maTimnm  change  in  p  is  less  than  10“^  times  the  maximum  value  of  p  .  The 
convergence  rate  depends  on  the  terrain  slopes,  since  N(p)  is  generally  proportional  to 
the  square  of  the  slope.  This  limits  the  maximum  slope  to  about  unity,  although  a 
practical  limit  is  0.75  ,  since  the  convergence  becomes  unacceptably  slow.  For  slopes  of 
0.5,  the  procedure  requires  about  5-7  iterations  for  convergence. 

The  accuracy  of  the  pressure  solution  is  improved  by  the  removal  of  the 
hydrostatic  pressure  field.  This  is  accomplished  in  the  transformed  coordinate  model  by 
averaging  the  potential  temperature  field  horizontally,  i.e.,  temperatures  are  interpolated 
onto  a  horizontal  level  before  averaging;  the  average  value  is  then  computed  on  the 
transformed  grid  and  subtracted  from  the  local  value  when  used  in  the  w-equation.  This 
procedure  is  applied  at  each  timestep  to  minimize  the  hydrostatic  contribution. 

Boundary  conditions 


'  Lateral  boundary  conditions  are  periodic  for  all  variables  and  therefore  need  no 
further  description.  The  domain  should  be  sufficiently  large  that  the  artificial  periodicity 
does  not  influence  the  results,  and  several  calculations  were  made  to  check  the 
sufficiency  of  the  domain  size.  The  upper  boundary,  z  =  2? ,  is  a  rigid,  stress-free  lid  but  a 
Rayleigh  damping  is  applied  to  the  vertical  velocity  over  the  top  6  grid  levels.  The 
damping  increases  toward  the  boundary  and  provides  an  effective  non-reflecting 
boundary  for  the  propagating  gravity  waves.  Thus 

m  =  %  =T23  =2/3  =  2^  =0  at  ^  =  D  (5.18) 

Surface  boundary  conditions  in  the  transformed  coordinate  system  require  careful 
consideration,  however,  since  we  are  specifically  interested  in  the  surface  forces.  The 
normal  velocity  component  is  easily  specified,  namely 

0  =  0  at  C  =  0  (5.19) 

The  momentum  fluxes  are  obtained  from  a  surface  layer  analysis  using  Monin-Obukhov 
similarity  to  relate  the  surface  stress  to  the  velocity  at  the  first  grid  point.  For  a  flat 
surface,  this  gives  T13  and  T23  at  z  =  0  from  the  relations 
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(5.20) 


where  Ml  and  vi  are  the  velocity  components  at  the  first  grid  point  above  the  surface,  at 
height,  zi .  The  roughness  length  is  zo .  «•  is  the  friction  velocity,  and  L  is  the  Monin- 

Obukhov  length,  defined  as 

^  _  U*Tq 

fcgHo 

Here,  k  is  von  Karman's  constant  with  a  value  of  0.4,  and  Ho  is  the  surface  heat  flux. 
The  function  ^  is  a  standard  representation  of  the  Monin-Obukhov  profiles  (Businger  et 

al.,  1971). 

The  surface  stress  derivation  has  been  described  explicitly  for  the  flat  surface, 
since  we  shall  extend  the  definition  to  the  terrain  case  by  assuming  the  same  surface  layer 
relations  in  the  local  tangent  plane.  This  is  not  strictly  accurate  since  the  buoyancy 
effects  are  controlled  by  the  vertical  gravity  direction  and  do  not  simply  rotate  with  the 
inclined  surface.  However,  sufficiently  close  to  the  surface  the  profile  is  logarithmic  so  if 
the  lowest  grid  level  is  close  enough  to  the  surface  the  errors  will  be  negligible.  In  the 
calculations  reported  below,  the  first  grid  level  is  less  than  5m  above  the  surface.  Clark 
(1977)  uses  a  simplified  treatment  of  the  lower  boundary  condition,  considering  only  the 
two  vertical  flux  components  for  the  horizontal  momentum,  but  points  out  the  need  for  a 
more  careful  analysis.  The  formulation  developed  below  is  a  more  complete 
representation  of  the  surface  layer,  and  is  consistent  with  the  subgrid  stress  model 
employed  in  the  LES. 

We  can  use  the  flat  surface  relations  to  determine  the  tangential  stress  on  the 
ground,  but  we  must  rotate  the  stress  tensor  into  the  Cartesian  frame  since  the  momentum 
equations  use  the  Cartesian  velocity  components.  This  procedure  requires  a  definition  of 
the  full  tensor  in  the  local  frame  defined  by  the  tangent  plane  and  the  direction  of  the 
tangential  velocity,  which  we  take  as 
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(5.22) 


fo  0  n 


I  =M* 


0  0  0 

,1  0  0, 


with  a*  obtained  from  the  tangential  velocity  at  the  lowest  grid  level.  The  neglect  of 
the  other  tensor  components  is  consistent  with  the  subgrid  stress  model,  since  the 
diagonal  stress  components  become  isotropic  near  the  wall.  The  frame  rotations  do  not 
change  the  isotropic  part  of  the  tensor  and  it  can  therefore  be  neglected.  In  a  real  flow, 
the  stress  tensor  is  not  isotropic  near  a  wall  but  the  description  of  these  effects  requires  a 
more  sophisticated  subgrid  closure  model.  We  expect  the  tangential  stress  to  provide  the 
dominant  momentum  transfer,  but  the  detailed  Reynolds  stress  behavior  near  the  wall 
should  be  a  topic  for  future  study. 


When  the  stress  tensor  is  rotated  into  the  Cartesian  frame,  we  obtain 

Tij  =ui(uinj  +  ujni] 

where  u  is  the  unit  vector  aligned  with  the  wind  (parallel  to  the  surface)  and  a.  is  the  unit 
normal  to  the  surface. 


The  consistency  and  accuracy  of  the  surface  boundary  condition  was  tested  by 
means  of  a  simple  flow  solution  for  a  plane  channel.  The  solution  involves  a  slight 
modification  to  the  terrain  transformation  described  above,  so  that  the  flow  in  the  two 
geometries  shown  in  Figure  5-1  could  be  calculated.  The  turbulence  scale,  Amax .  was  set 
equal  to  H  so  that  there  are  no  resolved  eddies  and  we  compute  a  standard  mixing  length 
solution  for  the  channel  profile.  The  flow  equations  involve  no  buoyancy  or  CorioUs 
terms  but  are  driven  by  a  streamwise  pressure  gradient,  giving  a  horizontally 
homogeneous  solution  with  no  velocity  component  normal  to  the  walls.  By  arranging  for 
the  normal  distance  between  the  two  plates  to  be  identical,  the  two  solutions  should  be 
related  via  a  simple  rotation. 

Profiles  of  the  mean  velocity  and  the  subgrid  Reynolds  stress  components  for  the 
two  solutions  are  shown  in  Figure  5-2.  The  predictions  from  the  rotated  channel  have 
been  rotated  into  the  frame  of  the  Cartesian  solution  and  show  very  good  agreement 
throughout  the  flow.  Correct  prediction  of  the  subgrid  kinetic  energy  close  to  the  wall 
required  careful  treatment  of  the  shear  production  term  at  the  lowest  grid  location  where 
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Figure  5-1. 


Schematic  geometry  of  the  horizontal  and  inclined  channel. 
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Figure  5-2. 


Normalized  mean  velocity  and  subgrid  Reynolds  stress  profiles  from  the 
two  geometries  in  Figure  5-1.  SoUd  line  is  the  horizontal  channel  result, 
symbols  indicate  inclined  channel. 


velocity  derivatives  at  the  surface  must  be  estimated.  The  estimate  was  obtained  from  an 
extrapolation  of  the  derivative  in  the  flow  to  maintain  consistency  with  the  log-layer 
profile.  In  summary,  consideration  of  the  rotated  Reynolds  stress  shows  that  all  the 
tensor  components  are  important  at  the  surface,  and  the  idealized  solution  for  a  rotated 
channel  demonstrates  that  the  model  boundary  conditions  and  finite-difference 
approximations  provide  a  consistent  representation  of  the  log-layer  solution  on  an 
inclined  surface,  with  a  smooth  match  to  the  surface  conditions. 

The  surface  heat  flux  is  simpler  to  deal  with,  since  we  only  need  consider 
rotations  of  a  vector.  We  have  chosen  to  apply  a  constant  heat  flux  condition  at  the 
surface,  and  interpret  the  flux  as  being  specified  for  a  umt  surface  area.  This  implies  a 
slightly  higher  heat  input  through  the  distorted  terrain  surface  due  to  the  area  increase.  In 
fact,  a  zero  slope  condition  is  applied  to  the  horizontal  temperature  fluxes  at  the  surface, 
and  the  entire  normal  flux  is  input  through  the  /fs  component  at  ^  =0.  Zero  slope 
conditions  are  also  applied  to  the  turbulent  kinetic  energy  equation,  although  this  is  not 
critical  since  the  surface  layer  is  dominated  by  the  balance  between  shear  production  and 
the  dissipation. 


S3  FLOW  PARAMETERS  AND  FORCE  BALANCE. 


The  topography  used  in  this  study  is  a  simple  sinusoidal  ridge,  i.e. 


(5.24) 


where  A,  is  the  wavelength  and  ho  is  the  maximum  height  The  maximum  slope  of  this 
terrain  is  nho/A,in  the  x-direction,  and  the  flow  is  statistically  homogeneous  in  the  y- 
direction.  We  restrict  attention  to  this  idealized  terrain  for  our  initial  studies,  and  hope  to 
develop  a  general  framework  of  understanding  that  can  be  extended  to  more  complex 
terrain. 


The  specification  for  the  domain  of  integration  is  a  horizontal  square  of  dimension 
L,  and  a  vertical  range  from  z=h(x,y)  up  to  z=D.  All  of  the  integrations  reported  here  use 
L=£)=4km  and  were  chosen  to  represent  a  mixed  layer  depth,  z,-,  of  roughly  1km.  The 
vertical  grid  extends  far  enough  for  the  upper  wave-damping  layer  to  have  no  significant 
effect  on  the  flow.  The  numerical  grid  is  uniform  in  the  horizontal,  with  48  points  in  each 
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direction,  but  non-uniform  in  the  vertical.  Most  of  the  integrations  used  61  points  in  the 
vertical,  with  10m  spacing  near  the  surface,  40m  in  the  mixed  layer,  and  an  expanding 

mesh  above  ^=1200m. 

The  flow  was  initialized  with  the  geostrophic  wind  and  a  uniform  vertical 
temperature  gradient  of  5°Ckra-i  above  C=800m;  a  constant  temperature  and  a  specified 
'boundary  layer'  horizontal  wind  were  prescribed  below  ^=800m.  The  initial  conditions 
were  chosen  to  provide  an  initial  well-mixed  layer  with  an  average  boundary  layer 
velocity  close  to  the  final  steady  state.  The  initial  conditions  are  not  critical,  but  we  need 
to  maintain  an  appropriate  resolution  during  the  evolution  of  the  flow,  so  we  require  the 
initial  mean  fields  to  be  reasonably  close  to  the  final  state.  The  flow  is  mamtained 
through  the  constant  geostrophic  wind  vector  (Ug  ,  Vg  )  in  (5.5).  The  turbulence  is 
initiatpH  by  a  small  random  perturbation  of  the  temperature  field,  and  the  surface  heat 
flux  is  linearly  increased  from  zero  to  its  specified  value  over  the  first  lOOOs  of  the 
integration.  Vertical  velocities  are  initialized  to  zero,  but  the  pressure  solution  ensures 
zero  divergence  after  one  time-step  so  mass  continuity  is  assured.  This  procedure  does 
induce  an  initial  transient  response  but  this  does  not  persist  long  in  the  presence  of  the 
convective  eddies.  The  convective  circulations  generally  require  about  3(X)0s  to  reach 
statistical  equilibrium  (for  the  parameters  considered  in  these  calculations),  but  the 
overall  boundary  layer  profiles  take  much  longer  to  achieve  the  balance  between  Coriolis 
and  friction  effects.  The  details  of  the  initialization  procedure  are  virtually  eradicated 
during  the  initial  growth  of  the  convective  eddies. 

One  of  the  most  important  effects  of  non-uniform  terrain  on  the  large-scale 
atmospheric  flow  is  the  surface  momentum  flux  balance  over  small-scale  features.  For 
the  periodic  conditions  of  the  large-eddy  simulation  domain,  integration  of  the 
momentum  conservation  equation  over  the  entire  volume  yields 

where  Fx  and  Fy  ate  the  average  surface  momentum  fluxes  of  x-  and  y-momentum  per 
unit  area  (and  per  unit  fluid  density  since  we  have  ignored  the  density  in  the  Boussinesq 


(5.25a) 

(5.25b) 
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equation).  The  surface  momentum  fluxes  are  composed  of  a  shear  stress  component  and 
a  pressure  component,  so  that  Fx=Tx  +  Px  ♦  and  Fy=Ty  , where 


'^x=j2isi'^nni  +  '^i3n3)dS 

(5.26a) 

Ty=Xl{ri2ni  +  i:23n3)dS 

Lj 

(5.26b) 

(5.26c) 

Here  n  is  the  unit  normal  at  the  surface,  S  where  4=0-  Momentum  fluxes  at  the  lateral 
boundaries  cancel  due  to  periodicity,  and  the  stress-free  rigid  lid  at  ^=2?  implies  zero  flux 
through  the  upper  boundary. 


If  we  define  a  boundary  depth,  Zi  then  we  can  also  define  a  mean  boundary  layer 


velocity  deficit 

“s-“* 

ZiL, 

and  similarly  for  vb. 


(5.27) 


The  velocity  (wb,Vb)  will  be  representative  of  the  bulk  mean  boundary  layer 
velocity  provided  that  most  of  the  deficit  occurs  below  z  =  z,-  +  /jq  /  2 .  The  displacement, 
hQi2,  is  the  average  displacement  of  the  terrain.  Given  the  sharp  inversion  cutoff  for  the 
turbulent  stresses,  the  perturbations  above  zi  are  small  and  are  only  maintained  by  gravity 
wave  transports.  In  the  steady  state  situation,  a  simple  geostrophic  balance  is  obtained, 

fi,(vB-yi)  =  T^+Px 

A(“s-“s)  =  -r).  *^'^*'’* 

This  balance  between  the  Coriolis  force  and  the  momentum  flux  is  well-known  in 
oceanographic  transport,  e.g.  Gill  (1982).  We  use  this  steady-state  relation  to  determine 
the  validity  of  our  numerical  solution,  continuing  the  calculation  until  the  balance  is 
achieved  within  a  specified  tolerance. 


Integrations  have  been  made  for  a  range  of  external  parameters,  including  surface 
roughness  (zo).  terrain  wavelength  (A),  terrain  height  (ho),  and  flow  stability  and 
direction.  This  large  number  of  independent  parameters  prevents  a  detailed  investigation 
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of  the  entire  parameter  space,  but  some  variation  of  each  has  been  considered.  The  flow 
stability  is  measured  by  the  ratio  w^jUg,  where 

=u^  + 

and 

wi  =-lrH<A 

w*  is  the  usual  convective  velocity  scale  (Deardorff,  1972)  and  Ho  is  the  surface  flux  of 
potential  temperature.  The  nominal  inversion  height,  z,.  for  all  the  runs  presented  below 

is  1000m. 


5.4  SURFACE  FORCE  RESULTS. 

The  series  of  runs  was  designed  to  provide  information  for  a  range  of  atmospheric 
flow  conditions.  We  examine  the  various  effects  of  terrain  geometry,  surface  roughness, 
atmospheric  stability,  and  wind  direction,  although  the  breadth  of  the  parameter  space 
prevents  exhaustive  definition.  The  major  features  to  be  presented  will  be  the  changes  in 
the  surface  force  balance  and  mean  boundary  layer  flow  induced  by  the  terrain.  The 
standard  flow  situation  will  be  a  5ms-i  geostrophic  wind  perpendicular  to  the  ridges,  with 
a  surface  heat  flux  of  0.03®Cms-i  and  a  surface  roughness  of  Im;  variations  from  these 
values  will  be  specified  where  appropriate. 

5.4.1  Terrain  Geometry  Effects. 

Since  we  are  considering  only  sinusoidal  ridges,  there  are  two  defining  parameters 
which  we  choose  as  the  horizontal  wavelength,  A,  and  the  maximum  slope,  s=iih(^X.  We 
consider  wavelengths  of  1km,  2km,  and  4km  relative  to  the  normal  inversion  height  of 
1km,  and  slopes  up  to  0.5.  We  first  examine  the  mean  flow  characteristics  in  the  quasi¬ 
steady  state,  where  the  average  Coriolis  force  is  nearly  balanced  by  the  surface  stresses. 

Figure  5-3  shows  the  mean  flow  over  ridges  with  maximum  slope  of  0.5  and 
/l=2km.  The  mean  streamlines  show  a  clear  reversed  flow  in  the  valley  and  acceleration 
over  the  ridge  crests,  as  expected,  but  the  separated  flow  regions  in  the  two  valleys  are 
different  The  difference  between  the  two  waves  indicates  the  reliability  of  the  averaging 
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procedure  and  shows  that  the  reversed  flow  region  contains  large  variability.  The 
streamlines  show  little  evidence  of  any  disturbance  in  the  overlying  atmosphere;  the  flow 
separation  effectively  allows  the  boundary  layer  to  accommodate  the  vertical 
displacements  without  significantly  disturbing  the  stable  region.  The  transverse  velocity 
component,  v,  shows  a  relatively  well-mixed  profile  with  Uttle  response  to  the  local 
terrain.  This  is  to  be  expected,  since  the  transverse  velocity  is  forced  by  the  ConoUs 
terms  which  act  on  a  slow  timescale  in  comparison  with  the  hill  transit  times.  The  v- 
component  therefore  represents  an  average  response  and  does  not  exhibit  much 
correlation  with  the  local  terrain  variations. 

An  example  of  the  instantaneous  flow  structure  that  contributed  to  the  average 
fields  in  Figure  5-3  is  given  in  Figure  5-4.  The  near-surface  streaklines,  obtained  by 
calculating  trajectories  in  the  instantaneous  («,v)-field  at  10m  above  the  surface,  show  the 
distinct  separation  region  between  the  two  ridges.  The  predominant  flow  near  the  surface 
is  in  the  positive  y-direction,  except  over  the  hUl  crests  where  the  w-component  is 
accelerated.  As  noted  in  the  mean  flow,  the  v-component  is  relatively  constant  over  the 
ridges  so  the  low  velocity  reversed  flow  in  the  x-direction  results  in  flow  along  the  valley 
axis.  The  separation  and  reattachment  lines  are  readily  visible  in  the  instantaneous  plot 
over  most  of  the  domain.  The  lines  are  meandered  by  the  large  eddies  and  become  very 
indistinct  in  places,  but  the  mean  flow  is  clearly  evident 

The  instantaneous  vertical  velocity  field  at  ^=500m,  is  also  shown  in  Figure  5-4, 
and  resembles  the  mid-level  convective  velocity  field  over  flat  terrain.  There  is  some 
suggestion  of  alignment  of  the  convective  updrafts  but  no  visible  effect  of  surface  terrain 
variations.  The  statistics  of  the  velocity  fluctuations  are  described  in  Section  5.6. 

The  surface  forces  determine  many  general  features  of  the  boundary  layer,  as 
discussed  previously,  and  we  now  examine  the  variation  of  the  forces  with  slope  and 
wavelength.  We  attempt  to  calculate  steady-state  forces,  although  there  is  always  a  slow 
growth  of  the  boundary  layer  due  to  the  surface  heat  input  The  achievement  of  steady 
conditions  requires  long  integration  times  since  the  balance  between  the  Coriolis  terms 
and  the  Reynolds  stress  gradients  is  only  reached  after  times  of  0(f  The  actual 
requirement  depends  on  the  choice  of  initial  conditions,  but  we  typically  integrate  for 
about  4/-1.  i.e.  40,000s.  Over  this  period,  the  growth  in  Zi  is  significant,  so  we 
artificially  reduce  Zi  during  the  integration  to  maintain  an  inversion  height  of  about 
1000m  relative  to  the  mean  terrain.  The  procedure  simply  involves  resetting  the 
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Figure  5-3. 


Mean  flow  over  2km  ridges  with  slope  0.5  and  zcplm.  (a)  mean 
streamlines,  contour  interval  is  200m^s“^  for  positive  values,  lOm^s  ^  for 
negative  values  (shown  dashed),  (b)  mean  transverse  velocity  component, 
contour  interval  of  0.4ms“^. 


X  (km) 


Figure  5-4.  Instantaneous  flow  field  over  2km  ridges  with  slope  0.5  and  zo^m. 

(a)  flow  streaklines  at  C=10m.  (b)  vertical  velocity  component  at  C=500m, 
contour  interval  of  O.bms'i,  dashed  contours  denote  negative  values. 
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temperature  field  to  the  initial,  undisturbed  condition  without  changing  the  velocity 
fields.  This  introduces  a  transient  response  as  the  flow  readjusts  to  the  new  temperature 
distribution,  but  allows  the  mean  boundary  layer  profile  to  achieve  accurate  balance  over 
the  longer  timescale. 


Figure  5-5  illustrates  the  evolution  of  the  instantaneous  surface  forces  over  the 
duration  of  an  integration;  the  pressure  and  tangential  shear  stress  components  are  shown 
for  the  A=2km  terrain  with  a  slope  of  0.5  and  a  reduction  in  Zj  at  t=24000s.  The  forces 
are  averaged  over  the  domain  but  there  are  significant  temporal  fluctuations  on  the 
timescale  of  the  large  eddies.  The  surface  force  results  used  in  the  subsequent  analysis 
are  averaged  over  the  last  3000s  of  the  run.  There  is  a  visible  discontinuity  in  the  forces 


Figure  5-5.  Time  history  of  the  domain-averaged  surface  forces  for  flow  over  2km 
ridges  with  slope  0.5  and  zo=lni.  Forces  are  displayed  per  unit  horizontal 
surface  area  and  per  unit  fluid  density. 
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induced  by  resetting  the  temperature  field  at  f=24000s,  but  this  is  not  large  compared  with 
the  turbulent  fluctuations.  The  pressure  force  decreases  foUowing  the  reduction  in 
inversion  height,  but  the  forces  recover  and  reach  a  reasonably  steady  state.  The 
integrations  were  generally  continued  until  the  balanced  state  (5.28)  was  achieved  with  a 
tolerance  of  10%,  i.e.,  the  difference  between  the  left  and  right  sides  of  (5.28)  was  less 

than  10%  of  the  total  force  magnitude  +Fyj  . 

The  variation  of  the  surface  forces  with  terrain  slope  for  the  2km  wavelength 
ridges  is  shown  in  Figure  5-6(a).  The  trends  are  clearly  evident;  the  pressure  force 
increases  with  slope,  as  does  the  lateral  stress,  Ty,  whtie  the  streamwise  surface  stress,  T^, 
tends  to  decrease  with  increasing  slope.  The  pressure  force  is  approximately  proportional 
to  the  square  of  the  slope;  each  doubling  of  the  slope  gives  roughly  a  four-fold  increase  m 
Px.  This  is  a  weU  known  result  for  small  slopes  in  the  linear  perturbation  regime,  but  is 
accurately  maintained  in  the  nonlinear  flow  with  a  slope  of  0.5.  However,  there  seems  to 
be  a  somewhat  fortuitous  balance  of  conflicting  trends  responsible  for  this  quadratic 
behavior  since  there  are  significant  changes  in  the  boundary  layer  flow  for  the  larger 
slopes.  The  variation  of  the  mean  boundary  layer  velocities,  (mb,  vb)  as  defined  in  (5.27), 
is  shown  in  Figure  5-6(b),  where  we  have  used  a  value  of  z/=950m  for  the  calculation. 
The  large  increase  in  surface  drag  over  steep  terrain  has  reduced  the  mean  boundary  layer 
velocity  from  4.2  ms'^  with  no  topography  to  3.3  ms'i  in  the  x-direction,  while  the  v- 
component  increases  from  1.8  ms-i  to  3.5  ms'i.  This  increased  deflection  of  the 
boundary  layer  flow  due  to  the  CorioUs  effect  is  responsible  for  the  increase  in  Ty  seen  m 

Figure  5-6(a). 

The  proportionaUty  between  Px  and  the  square  of  the  slope  is  surprising  in  light  of 
the  reduction  in  mb  ,  since  we  would  expect  the  pressure  to  vary  with  the  square  of  the 
streamwise  velocity.  However,  the  increasing  height  of  the  terrain  with  slope  mtroduces 
two  flow  effects  which  act  to  increase  the  pressure  force.  First,  there  is  some  tendency 
for  the  mean  velocity  to  increase  with  height  above  the  surface,  so  that  a  higher  terrain 
will  be  exposed  to  a  higher  velocity  for  a  fixed  bulk  mean  velocity.  Second,  the  capping 
inversion  acts  almost  like  a  rigid  lid,  as  shown  in  the  mean  flow  in  Figure  5-3.  This 
impUes  an  increased  speed-up  for  higher  terrain  due  to  the  constriction  effect,  and  a 
correspondingly  larger  pressiuB  force.  The  combination  of  these  conflicting  effects  to 
produce  a  quadratic  variation  with  slope  must  be  regarded  as  somewhat  fortuitous.  We 
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Figure  5-6.  Variation  of  surface  forces  and  mean  boundary  layer  velocity  with  terrain 
slope  for  2km  ridges  and  zo=lm.  (a)  average  pressure  force  (Px)>  ^nd 
tangential  stress  forces  {Tx ,  Ty)  per  unit  area,  (b)  mean  boundary  layer 
velocity  components  (us,  vg). 
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shall  return  to  the  problem  of  parameterizing  the  forces  on  ridges  after  we  have  examined 
the  results  for  other  parameter  variations. 

The  flow  dependence  on  the  height  of  the  terrain,  measured  by  the  ratio,  ho/zi,  is 
more  clearly  examined  by  varying  the  wavelength  of  the  ridges  but  holding  the  slope 
constant  Runs  were  made  with  wavelengths  of  1, 2,  and  4km  at  a  maximum  slope  of  0.5 
and  surface  roughness  length  of  Im.  This  implies  maximum  terrain  heights,  ho  of  159m, 
318m,  and  636m  respectively.  The  mean  streamlines  for  the  4km  wavelength  ridge  are 
shown  in  Figure  5-7.  It  is  clear  that  the  blockage  effect  of  the  ridge  induces  a  significant 
acceleration  of  the  flow  over  the  crest  relative  to  that  over  the  valley,  although  the  large 
separation  zone  reduces  the  effective  vertical  displacement  of  the  streamlines.  The  tu- 
component  of  the  surface  stress  is  shown  in  Figure  5-8  for  all  three  wavelengths,  and 
illustrates  the  variation  of  the  stress  perturbation  with  wavelength.  The  separation  region 
is  more  extensive  for  the  4km  ridges  but  the  peak  stress  at  the  crest  is  actually  lower  than 
the  shorter  wavelengths.  This  is  a  result  of  the  reduction  in  ub  in  combination  with  a 
larger  constriction  effect  for  the  636m  hill.  Profiles  of  the  mean  u-component  over  the 


Figure  5-7.  Mean  streamlines  for  flow  over  a  4km  ridge  with  slope  0.5  and  zo-lni. 
Contour  interval  is  lOOm^s"^. 
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Figure  5-8.  Surface  variation  of  the  T13  stress  component  for  three  terrain  wavelengths 

with  slope  0.5  and  zo=ltn..  ”  — —  X  =  1km, 

- A  =  2km, - A  =  4km. 

hill  crest  for  the  three  wavelengths  are  shown  in  Figure  5-9.  The  shorter  wavelength  hills 
produce  an  accelerated  flow  close  to  the  surface,  which  is  responsible  for  the  larger 
surface  stress  at  the  crest,  while  the  4km  wavelength  shows  no  evidence  of  a  low-level 
maximum  but  has  a  higher  speed  in  the  main  boundary  layer.  The  limited  vertical  extent 
of  the  aerodynamic  pressure  gradient,  being  controlled  by  the  horizontal  wavelength,  is 
probably  responsible  for  the  localized  acceleration  over  the  shorter  wavelength  ridges. 

The  variation  of  the  area-averaged  surface  forces  is  illustrated  in  Figure  5-10.  The 
pressure  force  increases  significantly  with  wavelength  as  a  result  of  the  increased  velocity 
over  the  crests  of  the  higher  ridges,  while  the  shear  stress  force  decreases  due  to  the 
enhanced  sheltering  effect  The  total  force  in  the  x-direction  increases  with  wavelength. 
These  results  suggest  a  parameterization  of  the  forces  in  terms  of  effective  botmdary  layer 
velocities,  and  we  shall  pursue  the  specification  of  the  appropriate  velocity  in  Section  5.5. 
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Figure  5-9.  Vertical  profiles  of  the  mean  velocity  component  normal  to  the  ridge  at  the 
crest  for  three  terrain  wavelengths  with  slope  0.5  and  zo=lni-  The  vertical 
scale  is  normalized  by  the  inversion  depth  to  facilitate  comparison. 

- A  =  lkm, - A  =  2km, 

- A  =  4km. 


5.4.2  Variation  With  Surface  Roughness. 

The  character  of  the  local  siuface  affects  the  details  of  the  surface  layer  profile  and  has  a 
strong  influence  on  the  surface  momentum  flux  for  flat  terrain.  In  the  case  of  free 
convection,  however,  surface  roughness  has  a  limited  influence  in  the  near-surface 
region;  the  bulk  of  the  boundary  layer  turbulence  is  unaffected  (Schinnarm,  1988;  Sykes, 
Henn  and  Lewellen,  1993).  We  examine  the  flow  over  2km  wavelength  ridges  at  a  slope 
of  0.5  with  zo=0.01m,  0.1m,  and  Im,  to  determine  the  response  of  the  flow  field  and  the 
surface  forces. 
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Figure  5-10.  Variation  of  average  surface  forces  with  terrain  wavelength  for  slope  0.5 
and  zo=lni- 


Figure  5-11  presents  the  streamwise  surface  stress  component  normalized  by  the 
flat  terrain  value.  To* ,  and  shows  the  effects  of  surface  roughness.  Tqx  is  the  value  of  T, 
for  flow  over  flat  terrain  with  the  same  geostrophic  wind  and  surface  heat  flux.  Flow 
separation  is  sensitive  to  surface  roughness,  with  a  pronounced  region  of  reversed  flow 
over  the  zo=lm  surface  but  no  separation  over  the  lower  roughness  terrain.  The  peak 
stresses  are  also  reduced  over  the  large  roughness  ridges,  but  this  is  probably  due  to  the 
reduction  in  ub  in  this  case.  The  trend  with  roughness  seems  consistent  with  the 
dependence  derived  by  Hunt  et  al.  (1988)  for  neutral  flow,  where  shear  effects  were 
shown  to  enhance  the  surface  stress  response  for  larger  roughness  lengths.  A  slope  of  0.5 
is  relatively  steep,  and  the  attached  flow  for  the  smaller  roughness  cases  is  probably  due 
to  the  mixing  effect  of  the  unstable  temperature  profiles.  The  turbulent  eddy  viscosity  in 
the  convective  layer  is  larger  than  the  neutral  case,  and  therefore  reduces  the  effective 
Reynolds  number  of  the  flow  and  inhibits  separation. 
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Figure  5-11. 


Surface  variation  of  the  T13  stress  component  for  three  surface  roughness 

lengths  and  2km  ridges  with  slope  O.5.. - zo  =  0.01m, 

- zo  =  0.1m, - zo  =  1®. 


The  surface  forces,  normalized  by  the  equivalent  flat  terrain  force,  are  shown  in 
Figure  5-12.  The  normalized  pressure  force  is  strikingly  independent  of  zo  for  slopes  of 
both  0.25  and  0.5,  while  the  shear  stress  force  shows  a  slight  reduction  at  the  higher 
roughness  value.  This  result  is  in  contrast  to  the  various  model  studies  for  neutral  flow 
over  ridges  (see  e.g.,  Taylor  et  al.,  1989),  which  show  a  significant  dependence  on 
roughness.  The  roughness  effect,  however,  is  thought  to  be  due  to  shear-enhancement 
(Hunt  et  al.,  1988)  and  is  not  present  in  the  asymptotic  theory  of  Sykes  (1980).  We 
hypothesize  that  the  shear  is  reduced  under  the  strongly  convective  conditions  of  the 
LES,  and  hence  there  is  less  variation  of  the  forces  with  surface  roughness.  We  should 
note,  however,  that  the  surface  stress  variation  shown  in  Figure  5-11  does  exhibit  a 
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Figure  5-12.  Normalized  surface  force  variation  with  roughness  for  flow  oyer  2km 
ridges  with  slope  0.5.  (a)  tangential  stress  component  in  the  x-direction. 
(b)  pressure  force.  SoUd  Une  and  symbol  is  slope  0.5,  dashed  line  and 
open  symbol  is  slope  0.25. 
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roughness  dependence.  We  have  performed  a  smaU  number  of  calculations  with  smaller 
surface  heat  flux  and  these  results  will  be  examined  in  the  next  section  to  determine  the 
effect  of  stability. 

5.4.3  Effect  of  Boundary  Layer  Stability. 

The  flow  parameters  of  the  previous  sections  are  generally  close  to  free 
convection  conditions.  The  average  friction  velocity  for  the  flat  terrain  cases  was 
between  0.28ms-i  and  0.4ms-i  while  the  convective  velocity  scale,  w*=lms“^  A  number 
of  integrations  were  made  with  w*=0.5ms“^  to  examine  the  dependence  on  stability. 
These  flows  still  have  significant  buoyancy  effects,  so  are  expected  to  be  different  from 
the  neutral  case,  but  we  expect  shear-induced  phenomena  to  be  comparable  with  the 
convective  dynamics. 

A  series  of  runs  was  made  with  A  =  1km  and  a  maximum  slope  of  0.25  to 
investigate  the  effects  of  surface  roughness.  The  surface  force  results  are  presented  in 
Figure  5-13,  normalized  by  the  flat  surface  x-stress  value,  and  show  close  agreement  with 
the  higher  heat  flux  force  values.  The  flat  stress  value  is  0.119m2s-2,  as  compared  with 
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Figure  5-14.  Surface  force  variation  with  wavelength  for  flow  over  ridges  with  small 
surface  heat  flux,  Ho=0.0045*Cms-i,  and  zo=lm.  SoUd  line  and  symbol  is 
Tx  ITqx,  dashed  line  and  open  symbol  is  P*  /Tqx,  both  for  slope  0.5.  Open 
squa^  are  Px  /Tqx  for  slope  0.25;  open  square  with  cross  is 

(Px  +UW  waveyTojc  for  slope  0.25. 

0.176m2s-2  for  Ho=0.03’Cms-K  but  the  pressure  force  shows  the  same  normalized  value 
and  is  also  relatively  constant  with  surface  roughness.  Thus,  with  an  unstable  layer 
characterized  by  L  approximately  equal  to  Zi,  we  obtain  the  same  qualitative  result  as  with 
L  much  less  than  z,-.  The  expected  dependence  on  roughness  for  neutral  flow  is  evidently 
not  manifested  until  the  stability  is  further  reduced.  It  is  possible  that  the  LES  may  not  be 
able  to  resolve  the  neutral  surface  layer  adequately  (Mason  and  Thompson,  1992)  and 
that  our  low  heat  flux  results  are  affected  by  this  deficiency.  However,  our  LES  results 
are  self-consistent  and  are  certainly  modified  by  the  buoyancy  effects,  preventing  direct 
application  of  the  neutral  layer  predictions.  A  detailed  simulation  of  the  neutral  layer  and 
its  transition  to  a  buoyancy  dominated  flow  is  beyond  the  scope  of  the  present 
investigation,  but  will  be  necessary  to  complete  the  description  of  stability  effects. 

The  variation  of  the  surface  forces  with  terrain  wavelength  and  slope  is  shown  in 
Figure  5-14,  and  indicates  an  anomalous  pressure  force  at  A  =  4km  and  slope  0.25.  The 
forces  are  normalized  by  the  flat  terrain  value  and  generally  compare  very  well  with  the 
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higher  heat  flux  results,  but  the  normalized  pressure  force  at  A  =  4km  and  a  slope  of  0.25 
is  almost  twice  the  expected  value.  The  reason  for  this  anomaly  is  immediately  clear 
upon  examination  of  the  mean  flow,  as  presented  in  Figure  5-15.  The  boundary  layer 
features  are  very  similar  to  previously  described  results,  but  there  is  a  very  distinct  steady 
gravity  wave  above  the  inversion.  The  vertical  profile  of  horizontally-averaged 
momentum  flux,  uw  ,  is  shown  in  Figure  5-16  and  indicates  a  gravity  wave  stress  of 
0.043m2s“2  due  to  the  standing  wave;  this  can  be  compared  with  the  total  surface  pressure 
force  stress  of  0.  lObmV^.  Contributions  from  the  other  components,  that  is,  the  resolved 
and  subgrid  fluctuation  correlations,  ate  seen  to  be  negligible  above  the  inversion.  When 
the  wave  stress  value  is  subtracted  from  the  total  pressure  force  for  the  A  =  4km  case  the 
normalized  force  matches  very  closely  with  the  higher  heat  flux  cases;  the  boundary  layer 
pressure  force  contribution  is  indicated  by  the  cross  in  the  open  square  in  Figure  5-14. 
This  is  the  appropriate  force  for  determining  the  mean  boundary  layer  velocity  profiles, 
because  the  relationship  (5.28)  implicitly  assumes  zero  stress  at  the  top  of  the  turbulent 
layer.  A  more  general  integral  of  the  momentum  equation  gives 


where  F,-  represents  the  average  Reynolds  stress  at  the  inversion.  This  is  precisely  the 


gravity  wave  momentum  flux  above  the  mixed  layer,  so  the  difference  between  this  and 
the  surface  force  gives  the  correct  momentum  balance. 


It  is  interesting  that  only  one  case  generates  a  significant  gravity  wave  stress  in 
comparison  with  the  surface  pressure  force.  The  Brunt-Vaisala frequency,  A,  where 


is  a  cutoff  frequency  for  the  waves;  higher  frequency  forcing  of  the  stable  layer  produces 
evanescent  disturbances.  Thus,  the  critical  wavelength  for  a  wind  speed  of  5ms-i  and  a 
potential  temperature  gradient  of  5°Ckm"^,  InUIN  ,  is  approximately  2.5km.  It  is  not 
surprising,  therefore,  that  we  only  see  a  significant  wave  response  for  A  =  4km  but  clearly 
the  wavelength  is  not  the  only  factor  of  importance.  The  larger  amplitude  terrain, 
maximum  slope  of  0.5,  produces  a  strong  separation  between  the  ridges  and  therefore 
reduces  the  streamline  displacement  in  the  main  boundary  layer,  similar  to  that  in  Figure 
5-7.  The  separation  evidently  modifies  the  boundary  layer  flow  enough  to  cause  a  net 
reduction  in  the  inversion  displacement,  and  the  wave  stress  is  reduced  to  0.015m2s  ^  in 
comparison  with  a  surface  pressure  stress  of  0.2m2s”2.  It  is  also  noteworthy  that  the 
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Figure  5-15. 


Mean  flow  over  4km  ridges  with  slope  0.5  and  zo=l®  and  surface  heal 
flux, /fo=0.0045*Cms-i.  (a)  mean  x-component  of  velocity,  contour 
interval  of  0.5ms-i.  (b)  mean  vertical  velocity  component,  contour 
interval  of  O.lmsr^  and  potential  temperamre  field,  contour  interval  of 
0.5*K  (shown  as  dotted  lines). 


momentum  flux  (m^s"^) 


Figure  5-16.  Vertical  profiles  of  horizontally-average  momentum  flux  for  the  4km  ridge 
with  slope  0.25  and  low  surface  heat  flux.  Short  dashes  represent  the 
subgrid  component,  long  dash-short  dash  represents  the  conrtbution  from 
resolved  scale  fluctuations.  Solid  line  is  the  total  flux,  including  the 
contribution  from  the  steady  mean  flow. 

higher  heat  flux  calculations,  //()=0.03'’Cms~K  showed  much  smaller  wave  amplitudes  for 
both  slope  values.  The  wave  stresses  were  0.02m^s~^  and  O.Olm^s"^  for  the  s=0.25  and 
5=0.5  cases  respectively.  The  ability  of  the  convective  boundary  layer  to  absorb  the 
terrain  streamline  deflections  and  prevent  disturbance  of  the  inversion  is  clearly  an 
important  factor  determining  the  stable  layer  response. 

Results  for  a  more  unstable  boundary  layer  were  obtained  with  M^=2ms  ^  zo=lni» 
and  Ho=0.03'’Cms-i  .  Surface  forces  for  a  slope  of  0.5  are  shown  in  Figure  5-17  as  a 
function  of  terrain  wavelength;  the  forces  are  normalized  by  the  flat  surface  stress  value 
of  0.061m2s-2.  The  effective  drag  coefficient  for  the  flat  surface  is  significantly 
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increased  by  the  instability,  but  there  is  also  an  increase  in  the  normalized  pressure  force 
relative  to  the  Ug=5ms-^  cases  at  the  two  longer  wavelengths  (see  Figure  5-10).  There  is 
evidence  of  the  increased  importance  of  buoyancy  in  the  A  =  4km  flows,  as  can  be  seen  in 
the  mean  streamlines  in  Figure  5-18.  The  separation  region  is  clearly  enhanced  by  the 
buoyant  instability,  accelerating  the  reversed  flow  up  the  lee  slope  of  the  ridge  and 
building  an  extensive  recirculation  region  over  the  valley.  These  changes  in  effective 
boimdary  layer  stability  are  enhanced  by  the  reduction  in  ub  due  to  the  surface  drag,  and 
are  therefore  more  important  for  the  longer  wavelengths. 

5.4.4  Variation  With  Flow  Direction. 

The  previous  sections  have  all  dealt  with  a  geostrophic  wind  normal  to  the  ridge 
axis.  There  is  some  interaction  between  the  two  horizontal  velocity  components  through 
the  Coriolis  effects,  so  we  examine  the  role  of  the  geostrophic  wind  direction  in  this 
section.  Four  wind  directions  were  calculated  for  A  =  2km  ,  M^=5ms  ^,zo— Itn, 
/fo=0.03*Cms-i  ^  and  a  slope  of  0.5.  A  direction  of  0“  is  the  standard  geostrophic  flow 


Figure  5-17.  Normalized  surface  force  variation  with  terrain  wavelength  for  flow  over 
ridges  with  slope  0.5  and  zo=lni,  at  low  geostrophic  wind  speed  of  2ms-i. 
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normal  to  the  ridges,  ±45*  implies  a  geostrophic  flow  of  ( 3.5msr^  ±3.5ms“^),  while  90 
represents  a  flow  of  (0, 5msri)  ^ 

The  variation  of  the  forces  with  direction  is  shown  in  Figure  5- 19(a).  The  force 
variation  seems  complicated  when  viewed  as  a  function  of  geostrophic  direction,  but 
shows  some  correspondence  to  the  mean  boundary  layer  velocity  (mb,  vb),  which  is 
displayed  in  Figure  5- 19(b).  The  boundary  layer  velocity  basically  responds  to  the 
rotation  of  the  geostrophic  wind,  but  shows  a  phase  lag  corresponding  to  the  drag- 
induced  deflection.  Although  the  general  behavior  is  as  anticipated,  with  an  approximate 
45°  turning  of  the  wind  in  the  boundary  layer  and  a  corresponding  response  in  the  surface 
forces,  the  quantitative  details  are  not  predictable  by  a  simple  rotation.  First,  there  is  a 
clear  anisotropy  in  direction  with  x-forces  behaving  differently  from  y-forces  since  there 
are  no  terrain  slopes  in  the  y-direction.  Also,  the  relation  between  the  pressure  force,  Px, 
and  the  x-component  of  velocity,  mb,  is  not  straightforward.  The  four  directions  chosen 
for  this  study  produce  two  pairs  of  mb  values.  For  simUar  values  of  mb,  the  pressure  force 
is  similar,  even  though  vb  is  significantly  different  However,  the  lower  value  of  mb 
produces  a  larger  pressure  force  than  expected  from  a  quadratic  law.  It  is  likely  that  this 
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behavior  is  due  to  the  enhanced  instability  at  low  values  of  mb  as  seen  in  the  previous 
section,  which  apparently  transfers  momentum  more  efficiently. 

Mean  velocity  fields  for  the  two  cases  with  geostrophic  wind  at  an  angle  of  ±45 
relative  to  the  ridge,  that  is,  geostrophic  winds  of  (  3.5ms“^  ±3.5ms"0.  are  shown  in 
Figure  5-20.  These  flows  are  quite  different,  with  mb  being  very  small  in  the  first  and  vb 
small  in  the  second.  However,  it  is  interesting  that  both  cases  show  that  there  is 
significant  wind  shear  in  the  component  with  the  small  layer-averaged  value.  The  major 
component  exhibits  a  well-mixed  region  above  the  hill  crests,  but  there  is  roughly  2ms  ^ 
shear  across  the  boundary  layer  depth  transverse  to  the  mean  flow. 

The  pressure  force  for  the  1ow-mb  values,  that  is,  geostrophic  wind  directions  of 
90°  and  45°,  is  much  larger  than  would  be  expected  by  scaling  the  0°  result  by  the  square 
of  Mb  .  However,  the  value  is  very  similar  to  that  obtained  with  Mg=2ms“^  and  a  similar 
value  for  mb  .  The  pressure  forces  scale  reasonably  well  with  the  value  of  the  boundary 
layer  velocity  if  we  use  the  stability-dependent  stress  values  based  on  mb  •  as  determined 
from  the  previous  section.  Apparently,  the  pressure  force  depends  on  the  effective 
stability  of  the  flow  as  determined  by  the  actual  mean  boundary  layer  velocity  normal  to 
the  ridges  and  is  largely  independent  of  the  flow  parallel  to  the  crests.  The  force 
parameterizations  are  discussed  in  the  next  section. 


5.5  SURFACE  FORCE  PARAMETERIZATION. 

The  extensive  parameter  studies  described  in  the  previous  section  can  be  used  as  a 
basis  for  a  general  surface  drag  parameterization  scheme  for  convective  flow  over  ridges. 
The  objective  of  such  a  scheme  is  the  representation  of  the  average  surface  forces  in 
terms  of  large  scale  parameters,  such  as  the  wind  speed,  surface  heat  flux,  terrain 
geometry,  etc.  We  postulate  a  functional  dependence  for  the  forces,  and  use  the  LES 
results  to  provide  quantitative  input  for  the  parameterizations. 

5.5.1  Pressure  Force. 


We  first  consider  the  pressure  force,  since  this  is  the  major  contributor  to  the  force 
variability.  The  general  expectation  is  that  the  pressure  force  will  vary  in  proportion  to 
the  square  of  the  terrain  slope  and  also  like  the  square  of  some  appropriate  flow  speed. 


85 


Figure  5-19.  Variation  of  surface  forces  and  mean  boundary  layer  velocity  with 
geos  trophic  wind  direction  for  2km  ridges  and  z(plm.  (a)  average 
pressure  force  (PJ,  and  tangential  stress  forces  (7^  .  Ty)  per  unit  area, 
(b)  mean  boundary  layer  velocity  components  (mb,  vb). 


86 


The  slope  dependence  is  anticipated  from  linearized  analyses  and  is  also  supported  by  the 
results  of  section  5.4.1.  We  suggest  that  the  bulk  average  boundary  layer  velocity,  ub  ,  is 
a  plausible  velocity  scale,  since  this  represents  the  mean  boundary  layer  flow  normal  to 
the  ridge.  The  first  approximation  for  Px  can  therefore  be  written  as 


Px  =  Cl 


(5.30) 


where  s  is  the  maximum  slope,  and  Ci  is  a  dimensionless  coefficient  We  do  not  expect 
Cl  to  be  a  constant;  it  must  describe  effects  other  than  flow  speed  and  terrain  slope.  The 
functional  dependence  of  Ci  on  other  flow  parameters  must  now  be  determined. 


The  surface  roughness  dependence  shown  in  Section  5.4.2  indicated  that  the 
pressure  force  was  closely  proportional  to  the  equivalent  flat  surface  stress,  if  all  other 
parameters  remained  fixed.  Since  we  wish  to  use  the  effective  boundary  layer  speed,  ub  , 
as  the  velocity  scale,  we  need  to  define  a  flat  surface  drag  coefficient,  cdo  .  such  that 

(t^+T^)''^  =  cdq{uI  +  vI) 

for  the  flat  surface.  Here,  (us ,  vb  )  is  the  steady-state  average  boundary  layer  velocity  for 
the  flat  surface,  and  the  force  dependence  is  chosen  to  be  consistent  with  horizontal 
isotropy.  We  have  not  used  the  x-component  of  the  stress  in  this  definition  of  cdo,  so  we 
are  not  entirely  consistent  with  the  force  normalization  employed  in  the  previous  section. 
However,  the  additional  complication  involved  in  defining  generalized  drag  coefficients 
for  the  different  direction  components  does  not  seem  warranted  at  this  stage.  From 
(5.31),  CDO  will  be  a  function  of  zo,  geostrophic  wind  speed,  Ug  ,  and  the  boundary  layer 
stability,  measured  by  w*.  Values  for  cdo  determined  from  the  LES  calculations  are 
shown  in  Figure  5-21  as  a  function  of  the  stability  parameter, 

“  /  2  2\1/2 

(4+n) 

Our  representation  of  the  pressure  force  now  becomes 
Px  =  C2CdqS^uI 

where  the  undetermined  coefficient,  C2 ,  accounts  for  other  flow  effects. 


The  different  terrain  wavelengths  in  Section  5.4.1  produced  different  pressure 
force  results,  so  this  must  be  accounted  for  in  C2  .  It  was  suggested  that  the  'constriction' 
effect  induced  by  the  stable  inversion  layer  could  be  responsible  for  the  increased  forces 
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Figure  5-21.  Bulk  mean  drag  coefficient,  cdo  » for  flow  over  flat  terrain  as  determined 
from  the  LES  calculations.  Variation  is  shown  for  different  roughness 
lengths  as  a  function  of  boundary  layer  stability,  Sb~ 

on  the  longer  wavelength  ridges.  Conservation  of  mass  implies  that  the  speed  over  the 
ridges  crest  must  increase  if  the  layer  depth  is  diminished,  and  we  observed  a  relatively 
horizontal  inversion  in  the  LES  results.  If  we  use  this  accelerated  speed  as  the  scaling 
velocity,  then  C2  contains  a  factor  inversely  proportional  to  the  square  of  the  depth  of  the 
boundary  layer  over  the  crest,  i.e., 

p  QcpqAI  (5.33) 

(1-^axAif 

where  /imax  and  z,-  ate  defined  relative  to  the  average  terrain  elevation. 

The  relation  (5.33),  with  C3  =  5,  provides  a  reasonably  good  description  of  the 
pressure  forces  for  most  of  the  calculated  cases,  as  can  be  seen  in  Figure  5-22.  The 
poorly  predicted  forces  are  associated  with  the  light  wind  cases  for  the  strongly  unstable 
boundary  layer,  and  for  two  of  the  wind  direction  values  from  Section  5.4.4.  Variations 
with  slope,  wavelength,  and  surface  roughness  are  well  predicted  for  the  moderately 
unstable  layer,  w*/Ug  =  0.2  ,  and  also  for  the  lower  heat  flux  cases,  where  w*/Ug  =  0.1 . 

It  was  postulated  that  the  anomalously  large  pressure  forces  shown  in  Figure  5-22 
were  due  to  the  effects  of  buoyant  instability  at  low  wind  speed.  It  was  suggested  in  the 
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Figure  5-22.  Comparison  between  LES  results  for  Px  and  the  estimate  from  (5.33). 


previous  section  that  the  stability  dependence  of  the  pressure  force  seemed  to  be 
determined  by  the  component  of  velocity  normal  to  the  ridge.  The  final  form  for  the 
pressure  force  is  therefore 


Px-^  ^DO  '“’in 


Z/  1“bU  (1-^/z.) 


(5.34) 


where  the  sign  of  the  boundary  layer  wind  is  now  accounted  for. 


The  functional  dependence  of  cpo  was  only  crudely  determined  from  the  LES 
results  shown  in  Figure  5-21,  and  simple  linear  interpolation/extrapolation  was  used  in 
the  appUcation  of  (5.34).  The  final  comparison  between  the  pressure  force  predicted  by 
(5.34)  and  the  calculated  LES  values  is  shown  in  Figure  5-23.  There  is  now  reasonably 
good  agreement  throughout  the  range  of  conditions  of  the  simulations,  with  only  one  case 
showing  an  error  of  more  than  about  20%.  This  is  the  largest  hill  at  the  smaUest  heat 
flux,  and  (5.34)  overpredicts  the  pressure  force  by  50%. 
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5.5^  Tangential  Stress  Forces. 

The  processes  controlling  the  tangential  turbulent  stresses  are  governed  by  the 
local  velocity  field  near  the  surface.  The  major  features  evident  in  the  force  variations  are 
a  decrease  in  Tx  due  to  the  flow  separation  between  the  ridges,  and  an  increase  in  Ty  due 
to  the  enhancement  in  the  transverse  velocity  component  An  obvious  approach  to  the 
parameterization  of  these  forces  is  therefore  to  assume  the  same  dependence  as  the  flat 
surface  case,  i.e., 

Tx  =  (4  +  4  f  ^  ^ 

=  (5.35b) 

and  we  only  need  to  determine  the  appropriate  drag  coefficient,  c/x)  •  Estimates  made 
using  the  same  value  as  in  the  pressure  force  (5.34)  are  compared  with  the  LES  results  in 
Figure  5-24.  There  is  a  clear  overprediction  of  the  tangential  forces  for  both  components, 
and  it  is  clear  that  these  components  do  not  show  the  stability  effects  observed  in  the 
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Ty  (LES)  (m2sr2) 

Figure  5-24.  Comparison  between  LES  results  for  (a)  Tx  and  (b)  Ty  and  the  estimate 
from  (5.35)  with  the  stability  dependent  coefficient  as  in  (5.34). 
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pressure  force  variations.  Revised  estimates  were  made  using  the  stabiUty  parameter  5b, 
based  on  the  mean  speed,  (ug  +  vlf^  .  rather  than  the  Ms-component,  and  the  results  ate 
shown  in  Figure  5-25.  This  also  contains  a  reduction  factor  on  7*  to  account  for  the 
separation  effects,  so  that  the  final  form  for  the  tangential  stress  is 


(5.36a) 

(5.36b) 


This  formulation  provides  a  reasonably  good  approximation  for  the  range  of  the  LES 
results. 


5.53  Extension  to  Complex  Terrain. 

The  results  of  the  previous  section  were  confined  to  sinusoidal  ridges,  but  the 
wide  rsmge  of  parameters  allowed  a  general  basis  for  the  force  parameterization.  In  order 
to  further  broaden  the  basis,  several  calculations  were  made  for  different  terrain 
geometry.  First,  variations  on  the  sinusoidal  shape  assumption  were  investigated  by 
extending  the  width  of  the  valley  region  in  comparison  with  the  ridge  itself.  The  terrain 
is  still  two-dimensional,  so  the  basic  parameterization  can  remain  as  described  by  (5.34). 
The  three  cases  are  illustrated  in  Figure  5-26,  showing  increasingly  isolated  ridges  but 
maintaining  the  Same  r.m.s.  terrain  slope.  The  ridges  become  higher  as  the  separation 
increases,  and  the  flow  perturbation  becomes  larger  but  is  localized  around  the  ridge.  The 
force  results  for  the  three  cases  are  given  in  Table  5-1. 


The  second  series  of  cases  are  three-dimensional  hills,  i.e.,  with  variation  in  the  y- 
direction.  A  modulation  of  the  ridge  terrain  was  tried,  with  elevation  defined  by  the 
relation 
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Thus,  h  =  ah2D  +  (1  “  oc)hiD  .where  h2D  is  the  ridge  topography  (5.24), 
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Ty  (LES)  (m2s-2) 

Figure  5-25.  Comparison  between  LES  results  for  (a)  Tx  and  (b)  Ty  and  the  estimate 
from  (5.36)  with  the  stability  dependent  coefficient  based  on  -i-  vg  j  . 


Figure  5-26.  Mean  w-velocity  for  convective  flow  over  four  ridge  shapes.  Standard 
sinusoidal  ridge  is  shown  at  top,  and  velocity  contours  are  in  ms"'. 
Geostrophic  wind  is  5nis~^ ,  the  surface  roughness  is  Im. 
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Table  5-1.  Forces  for  complex  terrain  shapes. 


Terrain 

ho 

(m) 

Max. 

slope 

Px 

(m^s"^) 

Py 

(m^s"^) 

Tx 

(m^s"^) 

Py 

3D-isolated 

318 

0.50 

0.025 

0.010 

0.158 

0.087 

0=1 

318 

0.50 

0.068 

0.036 

0.125 

0.094 

0=1 

159 

0.25 

0.018 

0.008 

0.164 

0.085 

Xjr*^\^21cni 

0=0.5 

318 

0.50 

0.121 

0.010 

0.1114 

0.130 

X^Xv=2km 

159 

0.25 

0.033 

0.002 

0.170 

0.092 

X,;(=2km 

0=0 

318 

0.50 

0.240 

0.0 

0.093 

0.162 

X,;(=2km,o=0 

(Fig.  5-26a) 

159 

0.25 

0.049 

0.0 

0.144 

0.098 

Fig.  (5-26b) 

225 

0.35 

0.057 

0.0 

0.138 

0.109 

Fig.  (5-26c) 

225 

0.35 

0.045 

0.0 

0.148 

0.099 

Fig.  (5-26d) 

318 

0.50 

0.094 

0.0 

0.148 

0.113 

and  a  is  an  interpolation  parameter  between  the  two  shapes.  Simulations  were  made  for 
the  parameters  given  in  Table  5-1,  with  values  of  a  of  0,  0.5,  and  1;  all  cases  used  a 
surface  roughness  height  of  Im  and  surface  heat  flux  of  0.03®Kms"^  The  table  also 
includes  an  isolated  circular  hill,  which  has  cosine  shape  with  a  radius  of  1km  in  the  4km 
square  domain.  Force  results  from  the  calculations  are  given  in  the  table. 

The  force  parameterization  given  in  (5.34)  is  clearly  a  specialization,  since  the 
only  component  of  the  pressure  force  is  in  the  x-direction  for  ridge  terrain.  We  require  a 
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general  directional  dependence  to  account  for  three-dimensional  terrain,  but  the  form 
must  collapse  to  (5.34)  for  two-dimensional  ridges.  The  generaliiation  of  the  r.m.s.  slope 
must  come  from  the  second-rank  tensor,  j,y,  defined  as 


_  dh  dh 

dxi  dxj 


The  force  parameterization  can  then  be  written  as 

SjjUBjUs 


Pi  =  10c£)o 

Ti  =cdq 


2/  “S 


Zi’UB 


(5.38a) 


(5.38b) 


where 

Ub  =  hflil 

had 

"  A(5) 

and  A(s)  is  the  maximum  eigenvalue  of  Sij.  This  somewhat  obscure  normalization  is 
required  to  ensure  that  the  force  is  proportional  to  the  modulus  of  the  boundary  layer 
velocity  for  horizontally  isotropic  terrain.  The  comparison  between  the  prediction  of  the 
pressure  force  from  (5.38)  and  the  complex  terrain  force  results  from  this  section  is 
shown  in  Figure  5-27.  The  results  are  encouraging,  and  suggest  that  the  given  form  may 
be  reasonably  accurate  for  a  wide  range  of  conditions. 


5.6  TURBULENCE  AND  DIFFUSION. 

The  presence  of  terrain  clearly  modifies  the  mean  velocity  field,  as  shown  in  the 
previous  section,  and  this  will  induce  changes  in  the  turbulence  fields.  The  turbulence  is 
distorted  by  the  additional  shears  and  responds  in  a  complicated  way.  Shear  generally 
produces  turbulent  kinetic  energy,  but  the  transient  shears  experienced  in  flow  over 
terrain  do  not  always  persist  long  enough  for  equilibrium  relationships  to  become 
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Figure  5-27.  Comparison  between  LES  pressure  forces  and  the  prediction  from  (5.38) 
for  the  range  of  terrain  shapes  described  in  Table  5-1  and  the  2d-ridges  in 

Figure  5-26. 


established.  The  dispersion  of  dust  particles  in  the  boundary  layer  will  be  affected  by 
both  the  mean  flow  and  turbulence  changes,  but  we  first  examine  the  turbulence  fields 
from  the  LES  calculations.  We  consider  the  sinusoidal  ridge  terrain  first 

5.6.1  Velocity  Variances. 

One  of  the  complicating  features  of  turbulent  flow  over  the  terrain  is  the  effect  of 
shear  distortion  on  the  turbulent  eddies.  The  stretching  and  tilting  of  the  eddies  as  they 
pass  over  the  terrain  is  a  non-equilibrium  process  if  the  timescales  of  the  distortion  are 
faster  than  the  eddy  turnover  timescales.  This  is  one  of  the  reasons  for  using  LES  to 
model  the  flow,  since  empirical  assumptions  about  the  Reynolds  stress  distributions  are 
unreliable  under  these  conditions.  We  do,  however,  use  an  equilibrium  assumption  for 
the  subgrid  stress  model,  which  assumes  that  the  timescale  of  the  subgrid  eddies  is  fast 
enough  that  they  remain  in  balance  with  the  local  wind  shear.  The  subgrid  turbulence 
scale,  Amaxj  is  generally  45m  for  the  calculations  reported  here,  so  that  the  subgrid 
turbulence  time  scale,  hJq,  is  roughly  45s  for  a  typical  r.m.s.  subgrid  turbulence  velocity, 
q,  of  lms”L  The  energy  dissipation  time  scale  is  longer  by  a  factor  of  which  is  4 
for  the  constants  used  in  the  turbulence  model  (Lewellen,  1977).  The  transit  time  for 
flow  over  the  ridges  depends  on  the  wavelength  and  mean  boundary  layer  velocity,  but 
most  of  the  results  are  for  a  2km  wavelength  and  a  geostrophic  wind  of  5ms“L  The 
average  boundary  layer  speed  is  somewhat  less  than  the  geostrophic  speed,  so  a  typical 
transit  timp.  for  passage  over  the  ridge  is  about  500s,  compared  with  the  dissipation  time 
of  about  200s.  The  subgrid  turbulence  consequently  responds  on  time  scales  faster  than 
the  distortion  time,  but  not  significantly  faster.  The  subgrid  perturbations  may  not  be 
completely  reliable,  therefore,  except  in  the  region  very  close  to  the  stuface  where  the 
scale  is  limited  by  distance  from  the  ground  and  response  time  scales  are  correspondingly 
faster.  The  subgrid  component  is  a  relatively  small  fraction  of  the  total  turbulence 
energy,  however,  so  inaccuracies  in  this  component  do  not  significantly  affect  the  overall 
results. 


The  general  behavior  of  all  three  components  of  the  velocity  variance  for  two 
different  terrain  slopes  is  presented  in  Figure  5-28.  Velocity  statistics  were  obtained  from 
a  transverse  average  across  the  integration  domain  and  a  3000s  time  average;  this  period 
was  determined  to  be  sufficiently  long  to  provide  reliable  statistics.  For  a  maximum 
terrain  slope  of  0.25,  the  variance  perturbations  are  relatively  small,  and  the  ambient 
vertical  profile  is  clearly  visible.  The  disturbance  near  the  surface  is  of  significant 
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amplitude  but  is  confined  to  a  layer  of  about  200m  depth.  The  amplitude  of  the 
disturbance  and  also  th^ertical  extent  is  increased  in  the  higher  slope  case.  The 
streamwise  component,  u'^ ,  shows  a  strong  increase  downstream  of  the  hill  crest  and 
above  the  surface.  This  is  a  region  of  high  vertical  shear,  as  the  surface  flow  is  reduced  in 
the  lee  and  is  associated  with  the  separation  bubble  boundary  for  the  slope  of  0.5.  The 
high  values  of  grow  from  a  very  shallow  layer,  originating  just  upstream  of  the  crest 
in  the  region  of  high  surface  stress  values,  and  spreading  vertically  downstream.  There  is 
evidence  of  a  reduction  in  at  some  height  above  the  surface  just  upstream  of  the  crest 
This  is  seen  as  a  lowering  of  the  contours  for  small  slope,  and  becomes  a  local  minimum 
at  a  slope  of  0.5.  The  vertical  velocity  variance,  w'^  ,  behaves  differently,  showing  an 
increase  above  the  surface  upstream  of  the  crest  but  smaller  values  in  th^shear  layer 
region  where  ^  is  ma^um.  The  most  dramatic  variations  are  seen  in  v'^ ,  however. 
There  is  an  increase  in  v'^  over  the  updope  of  the  ridge  that  becomes  very  large  for  the 
higher  slopes.  The  local  maximum  in  v'^  for  a  slope  of  0.5  is  almost  twice  as  large  as  the 
streamwise  fluctuation  variance  in  the  separating  shear  layer. 

The  variance  perturbations  are  qualitatively  consistent  with  existing  analyses  of 
turbulence  modification  in  flow  over  ridges.  Britter  et  al.  (1981),  discussing  the  response 
of  neutral  layer  turbulence,  suggest  that  the  variances  close  to  the  surface  are  governed  by 
the  tangential  stress  behavior,  which  maximizes  over  the  crest  Further  above  the  surface, 
the  turbulence  undergoes  rapid  distortion  effects  (Batchelor  and  Proudman,2954),  whmh 
reduces  u'^  over  the  crest  but  increases  the  two  transverse  components,  v'  and  w'  . 
Some  of  these  features  are  visible  in  the  convective  layer  results  in  Figure  5-28,  but  the 
existing  rapid  distortion  analyses  are  not  strictly  applicable  to  the  convective  layer 
situation.  The  analysis  has  been  carried  out  (Newley,  1985)  for  axisymmetric  turbulence 
with  ,  but  the  convective  layer  has  large  values  for  both  horizontal  components 

near  the  surface  and  a  small  vertical  velocity  variance  and  therefore  does  not  match  the 
assumptions  of  the  theory.  The  rapid  distortion  results  for  the  transverse  and  vertical 
components  are  sensitive  to  the  initial  turbulence  levels.  The  existing  analyses  also 
ignore  the  presence  of  the  solid  surface  in  their  derivation  of  the  velocity  spectrum  from 
the  vorticity  distortion.  The  more  complete  theory  of  Hunt  (1973)  would  be  required  for 
the  convective  situation  where  the  larger  eddies  ate  significantly  affected  by  the  wall.  It 
is  still  instructive  to  examine  the  quantitative  response  of  the  velocity  variances,  however, 
so  we  now  proceed  to  a  closer  examination  of  the  variance  profiles. 
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Vertical  profiles  of  the  three  variance  components  for  a  slope  of  0.125  are  shown 
in  Figure  5-29  in  comparison  with  the  flat  surface  result^  The  perturbations  are 
relatively  small,  but  several  features  are  clear.  First,  the  u'^  intensity  is  generally 
increased  throughout  the  layer,  indicative  of  a  nonlinear  response  due  to  the  repeated 
interaction  with  the  periodic  ridges  even  at  this  small  slope.  Relative  to  the  mean  value, 
u'^  is  increased  near  the  surface  but  is  reduced  from  the  flat  sxuface  value  over  the  crest 
of  the  ridge.  Above  the  valley,  is  increased  except  for  a  reduction  in  the  lowest 
100m.  The  transverse  component  shows  an  increase  in  the  lower  half  of  the  boundary 
layer  both  over  the  crest  and  the  valley,  while  the  vertical  component  shows  a  general 
reduction  except  for  the  lowest  100m.  The  results  must  be  viewed  with  caution  near 
the  surface,  since  the  bulk  of  the  variance  is  represented  by  the  subgrid  model  using  the 
'quasi-equilibrium'  assumption. 

The  trends  observed  in  the  small  slope  case  are  more  obvious  when  the  slope  is 
increased  to  0.5.  Figure  5-30  shows  the  variance  profiles  from  the  large  slope  calculation 
and  the  general  increase  in  all  three  components  is  immediately  obvious.  The  horizontal 
variance  components  are  roughly  doubled  in  the  main  boundary  layer,  with  significant 
perturbations  in  the  lower  altitudes,  u'^  shows  its  ma^um  amplitude  in  the  shear  layer 
that  separates  from  the  lee  slope  of  the  ridge,  while  shows  the  dramatic  increase  in 
the  lowest  200m,  as  noted  in  Figure  5-28.  The  variation  in  the  depth  of  the  boundary 
layer  relative  to  the  local  surface  elevation  is  clearly  visible  in  the  profiles,  with  a 
maximum  difference  of  about  300m  between  crest  and  valley. 

The  profiles  of  w'^  in  Figure  5-30  also  show  an  increase  in  the  main  boundary 
layer,  with  a  local  maximum  associated  with  the  separating  shear  layer.  There  is  a  close 
correspondence  between  w'^  and  u'^  in  the  relative  variation  in  the  main  boundary  layer, 
both  components  show  a  minimum  over  the  crest  and  a  maximum  over  the  valley.  The 
average  value  of  w'^  is  increased  much  less  than  the  horizontal  components,  however. 
This  is  consistent  with  the  general  observations  of  increased  horizontal  variance  in 
complex  terrain,  although  the  sinusoidal  ridges  demonstrate  that  the  topography  need  not 
be  complicated;  a  significant  slope  is  all  that  is  required  to  generate  the  large  vertical 
gradients  and  induce  high  amplitude  horizontal  velocity  fluctuations. 

A  more  striking  feature  of  the  vertical  velocity  variance  profiles  is  the  strong 
maximum  just  above  the  inversion.  The  elevated  maximum  is  all  resolved  scale  energy. 
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Figure  5-29.  Vertical  profiles  of  velocity  variances  for  max.  slope=0.125  at  the  hi 
crest  and  valley.  The  dashed  lines  are  for  the  corresponding  flat  case. 
Both  resolved  and  total  variances  are  shown. 
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Figure  5-30.  Vertical  profiles  of  velocity  variances  for  max.  slope=0.5  at  the  hill  crest 
and  vaUey.  The  dashed  lines  are  for  the  corresponding  flat  case.  Both 
resolved  and  total  variances  are  shown. 
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and  examination  of  the  instantaneous  flow  fields  shows  immediately  that  it  is  due  to  a 
wave  disturbance  in  the  stable  region.  The  profiles  in  Figure  5-30  are  obtained  from 
averages  in  the  y-direction  as  well  as  a  time  average,  and  a  typical  realization  of  the 
vertical  velocity  field  at  z=1200m  is  shown  in  Figure  5-31.  A  coherent  wave  pattern  is 
clear,  and  seems  to  be  associated  with  a  roll-like  convective  structure  in  the  boundary 
layer.  The  wave  is  aUgned  at  30"  to  the  left  of  the  geostrophic  wind  direction,  consistent 
with  the  mean  wind  in  the  boundary  layer.  The  Ekman  balance  between  the  mean 
boundary  layer  flow  and  the  surface  drag  forces  implies  that  large  slope  cases  produce 
large  angles  between  the  rolls  and  the  geostrophic  direction.  This  increases  the  shear  at 
the  inversion  and  produces  a  stronger  relative  flow  across  the  rolls,  so  that  smaller  slopes 
show  much  weaker  wave  disturbance.  The  trapped  wave  velocity  variance  is  an  unsteady 
feature  that  varies  with  the  inversion  structure  and  convection  pattern;  it  does  not  appear 
to  influence  the  turbulent  layer  significantly  so  the  dynamics  of  the  waves  will  not  be 
pursued  further. 

The  most  dramatic  feature  of  the  turbulence  response  within  the  boundary  layer  is 
undoubtedly  the  increase  in  near  the  surface.  We  are  not  aware  of  previous 
measurements  or  predictions  of  this  phenomenon,  and  we  therefore  examine  the 
mechanism  in  more  detail.  The  instantaneous  v-velocity  field  at  f=50m  is  shown  in 
Figure  5-32,  and  intense  regions  of  high  transverse  velocity  are  clearly  visible.  The 
regions  are  elongated  in  the  direction  of  the  near-surface  flow  and  are  located  on  the 
upslopes  of  the  ridges.  The  high  transverse  velocities  are  associated  with  high 
streamwise  vorticity,  as  can  be  seen  in  Figure  5-33.  The  x-component  of  vorticity  at 
^=80m  shows  very  intense  features  along  the  upwind  slopes,  while  the  vertical  section  at 
x=500m  (midway  up  the  slope)  shows  localized  vortices  above  the  surface  layer  of 
negative  vorticity.  The  surface  layer  is  produced  by  the  mean  v-component  of  velocity  in 
the  FVnian  layer,  but  the  isolated  vortices  are  clearly  visible  above  it  with  both  positive 
and  negative  signs  and  a  scale  of  about  250m. 


The  velocity  fluctuations  associated  with  these  intense  streamwise  vortices  are 
strongly  affected  by  the  presence  of  the  solid  boundary,  since  the  vortices  are  close  to  the 
wall.  This  is  the  reason  for  the  preferential  manifestation  of  the  kinetic  energy  in  the  v- 
component,  rather  than  equally  amongst  v  and  w.  An  idealized  reflective  boundary 
condition  at  the  wall  would  reduce  the  normal  component  of  the  velocity  to_^ro  and 
double  the  transverse  velocity  component,  giving  a  strong  enhancement  of  v'  at  the 
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Figure  5-32.  Transverse  velocity  at  a=50tn  for  max.  slopesO.5.  Contour  intervals  are 
0.5ms“i. 
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Figure  5-33. 


Streamwise  component  of  vorticity  for  max.  slope=0. 
areO.OOSs"^  (a)  at  ^=80m,  (b)  at  x=500m. 


5.  Contour 


surface.  Th£waU  layer  drag  reduces  the  velocity  very  close  to  the  surface,  but  the  large 
increase  in  on  the  lower  side  of  the  vortex  is  clearly  evident  in  the  LES  results. 

The  vortices  are  reminiscent  of  the  Taylor-GorUer  vortices  observed  in  boundary  layer 
flow  over  a  concave  surface,  e.g.,  Tani  (1962),  So  and  Mellor  (1975)  and  Hoffman  et.  al. 
(1985).  Although  the  main  boundary  layer  is  deep  compared  with  the  terrain  scales,  the 
velocity  shear  is  largely  confined  to  a  relatively  thin  surface  layer.  The  valley  region 
presents  a  concave  surface,  and  we  might  therefore  expect  production  of  coherent 
streamwise  vortices  in  the  local  flow.  However,  the  intensity  of  the  vortices  observed  m 
our  LES  calculations  appears  to  be  significandy  larger  than  experimental  observations  m 
curved  boundary  layer  flow.  This  is  due  to  enhancement  by  the  shear  distortion  in  the 
flow  over  the  ridges.  Experimental  flows  have  generally  involved  minimal  streamwise 
acceleration  in  an  attempt  to  isolate  the  curvature  mechanisms,  but  in  the  ndge  flow  the 
vorticity  production  occurs  in  a  region  of  strong  acceleration.  Streamwise  vorticity  is 
enhanced  by  stretching  of  the  vorticity  field  as  the  flow  accelerates  toward  the  crest  of  the 
ridge,  increasing  the  v-velocity  and  also  stabilizing  the  vortex.  The  rate  of  stretching  of 
vorticity  aligned  with  the  mean  flow  direction  is  given  by 

dxj 

and  the  spatial  distribution  of  this  quantity  is  illustrated  in  Figure  5-34. 
potential  flow  over  a  sinusoidal  surface,  the  stretching  rate  is  simply  2e  dh/dx 
where  h  is  the  surface  elevation  and  A  is  the  terrain  wavelength.  The  linearized  vertical 
damping  scale  is  therefore  about  300m  for  the  2km  ridges  and  the  maximum  stretching  is 
located  at  the  maximum  slope  on  the  upwind  side  of  the  ridge.  The  general  shape  of  the 
stretching  field  in  the  nonlinear  flow  over  the  steep  ridges  is  qualitatively  similar  to  the 
linearized  result,  and  the  stretching  rates  are  fast  enough  to  produce  significant  mcrease  m 

the  local  vorticity. 

The  persistence  of  the  vortices  is  illustrated  in  Figure  5-35,  which  shows  the  v- 
velocity  midway  up  the  ridge  slope  and  40m  above  the  surface  as  a  function  of  time  from 
the  two  runs  with  maximum  slopes  of  0.25  and  0.5.  The  sUce  across  the  4km  wide 
domain  shows  intermittent  structure  at  a  slope  0.25,  with  occasional  'streaks  of  high 
speed  flow  lasting  for  about  1000s.  At  a  slope  of  0.5,  however,  two  distinct  narrow 
regions  of  very  high  v-velocity  migrate  in  the  positive  y-direction  at  about  3ms-i  for  the 
whole  3000s  period.  The  magnitude  of  the  v-velocity  in  the  streak  is  modulated  over  the 
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Figure  5-34.  Stretching  rate  of  vorticity  aligned  with  the  mean  flow  for  max.  slope=0.5. 
Contour  increments  ate  5xl0~^s“i. 


time  period  but  the  identity  of  the  structure  is  maintained  throughout,  providing  strong 
evidence  of  the  stabUity  of  these  vortices.  The  difference  in  persistence  time  between  the 
two  slope  cases  seems  to  be  related  to  the  existence  of  a  reverse  flow  in  the  valley.  At  a 
slope  of  0.5,  the  separated  flow  region  is  much  mote  complex  than  the  classical  curved 
boundary  layer  of  the  Taylor-Gortler  instability,  and  streamwise  vorticity  is  probably 
generated  through  instabilities  of  the  separating  shear  layer  The  strong  stretching 
mechanism  is  centered  on  the  reattachment  point,  however,  and  vortices  are  therefore 
almost  fixed  relative  to  the  terrain  surface.  Once  initiated,  a  vortex  can  evidently  sustain 
itself  under  the  stretching  action  for  long  time  periods. 


The  vorticity  production  mechanism  discussed  above  is  controlled  by  the  Taylor- 
Gortler  instability  mechanism,  which  depends  on  the  velocity  profile  shape  over  the 
concave  valley  surface.  The  vertical  extent  of  the  instability  is  determined  by  the  profile 
shear  and  wiU  therefore  depend  on  surface  roughness.  The  turbulence  enhancement 
induced  by  the  vortex  stretching  mechanism,  on  the  other  hand,  scales  with  the 
wavelength  of  the  terrain  and  does  not  depend  strongly  on  ^tors  like  the  surface 
roughness  or  boundary  layer  stability.  Figure  5-36  shows  the  v'^  profiles  for  different 
roughness  lengths  at  the  mid-point  of  the  upslope.  It  can  be  seen  that  the  depth  of  the 
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Figure  5-35.  Transverse  velocity  at  x=500ra,  z=40m  as  a  function  of  time  and  y  for 
(a)  max.  slope=0.25,  (b)  max.  siope=0.5. 
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Figure  5-35.  Transverse  velocity  at  x=500ni,  z=40m  as  a  function  of  time  and  y  for 
(a)  max.  slope=0.25,  (b)  max.  slope=0.5  (Continued). 
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Figure  5-36.  Vertical  profiles  at  jc=500m  of  total  and  resolved  transverse  velocity 
variance  for  max.  slope=0.5.  SoUd  lines:  zo=lni;  long  dashes:  zo=0.1m; 
and  short  dashes:  2o=0.01m. 


enhanced  v'^  layer  does  decreases  with  roughness,  although  the  magnitude  also  appears 
to  be  roughness-dependent 

Figure  5-37  shows  the  turbulence  profiles  for  a  caseAvith  smaller  buoyancy  flux, 
but  other  parameters  as  in  Figure  5-30.  The  increase  in  near  the  surface  is  again 
obvious,  suggesting  the  importance  on  the  near-surface  shear  instabiUties  and  the 
stretching  effect,  since  neither  mechanism  depends  on  buoyancy  or  convective  eddies. 
The  flat  terrain  profiles  show  that  the  turbulence  is  not  dominated  by  convection  in  this 
case,  since  the  variances  maximize  close  to  the  ground. 
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Figure  5-37.  Vertical  profiles  at  ;c=500m  of  resolved  and  total  transverse  velocity 
variance  for  max.  slope=0.5,  /fo=0.0045°Cms“^  The  dashed  lines  are  for 
the  corresponding  flat  case. 

5.6.2  Heat  and  Momentum  Fluxes. 

The  off-diagonal  Reynolds  stresses  are  subject  to  the  same  rapid  distortion  effects 
as  the  velocity  variances,  so  that  the  eddy  viscosity  relation  is  only  valid  close  to  the 
surface.  The  vertical  profiles  of  momentum  and  heat  fluxes  from  the  small  slope  case 
with  7Vm  ridges  are  shown  in  Figure  5-38,  and  compared  with  the  flat  terrain  case.  The 
profiles  of  u'w'  and  v'w'  are  very  similar  in  character,  and  show  the  major  features  seen 
in  observational  studies.  The  Reynolds  stresses  shows  an  elevated  minimum  over  the 
crest  of  the  ridges,  with  a  maximum  at  the  surface  and  an  increase  from  the  flat  terrain 
profile  at  higher  altitudes.  The  elevated  u'w'  minimum  is  just  below  z=100m  while  the 
Vw'  minimum  is  around  200m.  Estimates  of  the  height  of  the  stress  minimum  from  the 
neutral  flow  analyses  of  Sykes  (1980)  and  Belcher  et  al.  (1993)  are  close  to  the  inner 
layer  scale  height,  and  the  theoretical  prediction  has  been  compared  with  field 
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Figure  5-38.  Vertical  profiles  of  vertical  fluxes  for  max.  slope=0.125  at  the  hill  crest 
and  valley.  The  dashed  lines  are  for  the  corresponding  flat  case.  Both 
resolved  and  total  fluxes  are  shown. 


observations  by  Belcher  et  al.  Using  the  definition  of  the  inner  scale  from  Belcher  et  al. 
and  assuming  the  horizontal  scale  of  the  sinusoidal  ridges  as  L=?JA,  the  inner  scale  is 
obtained  from  the  expression 

41n(4/zo)  =  2K-^^ 

as  /j  =  50m .  The  result  in  Figure  5-38  is  somewhat  higher  than  this  estimate  and  may  be 
affected  by  the  LES  resolution.  The  minimum  in  u'w'  is  located  at  the  point  where  the 
resolved  scale  component  begins  to  dominate,  and  could  possibly  move  downward  if  the 
resolution  was  increased.  The  agreement  with  the  theoretical  estimate  is  reasonably 
good,  however,  and  indicates  that  most  of  the  flux  transport  is  adequately  resolved  m  this 

calculation. 


The  momentum  flux  profiles  in  the  valley  show  the  opposite  phase  in  the 
perturbation  near  the  ground,  but  the  fluxes  are  increased  throughout  the  main  boundary 
layer.  The  fluxes  are  generally  higher  over  the  valley  at  heights  above  the  inner  layer 
scale,  and  are  similar  to  the  velocity  variances,  which  indicate  a  nonlinear  response  at  this 
slope  and  a  net  increase  in  the  turbulence  level.  The  increase  in  average  momentum  flux 
is  manifested  in  the  total  surface  drag,  discussed  in  Section  5.5,  which  shows  a  pressure 
force  of  roughly  10%  of  the  flat  stress  value. 

The  heat  flux  perturbations  show  a  simpler  response  than  the  Reynolds  stresses, 
with  a  reduction  above  the  crest  and  an  increase  over  the  valley.  The  surface  heat  flux  is 
a  prescribed  constant  for  the  flow,  so  there  can  be  no  net  change  over  the  ridges,  but  there 
is  also  no  indication  of  an  inner  layer  in  these  profiles. 

At  a  larger  slope,  the  flux  perturbations  are  much  more  significant  Figure  5-39 
shows  the  profiles  for  a  slope  of  0.5.  The  large  increase  in  surface  drag  is  reflected  m  the 
uV  profiles,  particularly  over  the  valley,  where  the  momentum  flux  across  the 
separating  shear  layer  is  clearly  the  dominant  transfer  mechanism.  The  profile  over  the 
crest  still  shows  the  elevated  minimum,  at  a  slightly  lower  elevation  than  the  small  slope 
case  and  a  similar  response  in  v'w'.  The  bulk  of  the  transfer  of  v-momentum  is  also 
located  near  the  shear  layer  above  the  valley,  but  the  flux  profiles  show  considerable 
oscillation  vsdth  altitude.  The  heat  flux  is  also  clearly  reduced  over  the  hill  crest  and 
augmented  above  the  valley.  The  full  two-dimensional  structure  of  the  fluxes  is 
illustrated  in  Figure  5-40,  which  shows  the  heat  flux  maximizing  on  the  lee  slope  of  the 
ridges;  apparently  the  convective  eddies  are  preferentially  initiated  in  the  upslope  flow  of 
the  separation  region.  The  heat  absorbed  on  the  windward  slope  is  not  sufficient  to 
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Vertical  profiles  of  vertical  fluxes  for  max.  slope=0.5  at  the  hill  crest  and 
valley.  The  dashed  lines  are  for  the  corresponding  flat  case.  Both 
resolved  and  total  fluxes  are  shown. 
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promote  upwind  separation,  but  produces  a  buoyant  separating  shear  layer  that  tends  to 
rise  and  form  a  convective  updraft,  v'w'  also  shows  the  strong  maximum  in  the  shear 
layer  that  continues  up  to  the  inversion,  producing  large  horizontal  variations  throughout 
the  depth  of  the  boundary  layer.  This  contrasts  with  u'w'  and  w'd',  which  become  more 
horizontally  uniform  above  the  hills. 

5.6  J  Dispersion  Effects. 

The  turbulence  levels  over  complex  terrain  are  generally  increased  relative  to  flat 
terrain.  Grant  and  Mason  (1989)  have  suggested  that  the  turbulence  above  the  crests  of 
the  hills  is  comparable  to  that  over  a  flat  surface  with  an  equivalent  surface  stress,  and 
there  is  rough  confirmation  of  this  hypothesis  in  the  LES  results.  As  we  have  seen,  the 
flow  field  is  nearly  homogeneous  in  the  horizontal  at  heights  more  than  about  100m 
above  the  ridge  crests;  this  statement  does  not  apply  to  flux  of  the  momentum  component 
along  the  ridge,  vV,  but  the  variance  profiles  and  u'w'  profiles  are  close  to 
homogeneous.  Turbulence  variances  do  not  show  exactly  the  same  shape  as  a  flat  terrain 
profile,  but  the  magnitudes  are  comparable.  We  should  therefore  expect  enhanced 
dispersion  over  complex  terrain,  and  the  appUcation  of  the  simple  modeling  scheme  of 
Section  2  should  produce  the  correct  trends.  The  total  surface  forces  obtained  from  (5.38) 
should  be  used  to  estimate  u*  in  the  parameterization.  However,  the  actual  dispersion 
process  is  more  complicated  than  a  simple  increase  in  turbulence  intensity. 


In  order  to  obtain  a  picture  of  the  dispersion  effects,  the  LES  calculation  was  used 
to  track  a  number  of  "massless"  Lagrangian  particles  for  a  period  of  time.  The  particles 
were  released  into  the  convective  boundary  layer  over  terrain  of  slope  0.5,  as  illustrated  in 
Figure  5-3.  This  is  a  separated  flow,  and  the  turbulence  intensities  are  shown  in  Figure 
5-28.  A  cluster  of  100  particles  were  released  into  the  fully  developed  turbulent  flow  at 
time  t=0,  say.  The  particles  were  initially  located  in  a  10x10  array  over  a  250m 
horizontal  square  at  an  altitude  of  400m.  The  square  was  positioned  above  the  crest  of 
the  ridge,  which  has  an  elevation  of  318m  for  this  case.  The  particles  are  transported  by 
the  resolved  scale  eddy  field,  and  no  subgrid  diffusion  effects  were  included. 

The  subsequent  positions  of  the  particles  at  15min  intervals  out  to  t=l  hour  are 
shown  in  Figure  5-41.  The  cluster  only  provides  a  rough  estimate  of  the  dispersion 
effects,  but  it  is  clear  that  there  is  a  rapid  spread  of  the  distribution.  The  mean  wind 
direction  is  clear  in  the  plan  view,  with  particles  moving  almost  45®  to  the  left  of  the 
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Figure  5-41.  Particle  dispersion  over  1  hour  from  a  localized  release  in  convective  flow 
over  periodic  ridges,  maximum  terrain  slope  is  0.5,  and  geostrophic  wind 
is  5ms“^  in  the  x-direction. 
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geostrophic  wind  direction,  i.e.,  the  positive  x-direction.  This  is  to  be  expected  from  the 
surface  force  results,  which  determine  the  mean  boundary  layer  velocity  through  the 
balance  with  the  Coriolis  force.  The  particles  spread  more  rapidly  in  the  x-direction,  and 
the  elevation  view  shows  that  the  spread  is  closely  associated  with  the  trapping  and 
retardation  of  particles  in  the  valleys.  The  flow  is  turbulent,  so  that  even  though  the  mean 
streamlines  are  closed  within  the  recirculation  region  there  is  a  strong  exchange  of  fluid 
between  the  valleys  and  the  overlying  boundary  layer.  The  spread  in  the  x-direction  is 
about  8km  over  the  1  hour  period,  while  that  in  the  y-direction  is  about  3km.  The  latter  is 
consistent  with  the  magnitude^  ,  as  shown  in  Figure  5-28,  but  the  x-spread  is  much 
larger  than  implied  by  the  distribution.  The  mean  flow  variations  are  clearly 
responsible  for  much  of  the  diffusion  in  the  x-direction,  and  this  must  be  accounted  for  in 
a  proper  description  of  dispersion  in  complex  terrain. 

The  mean  flow  effects  can  be  represented  in  late-time  dispersion  models  through 
appropriate  specification  of  the  mean  wind  profile  over  complex  terrain.  The  flow  field 
shown  in  Figure  5-3  can  be  averaged  in  the  horizontal,  and  will  yield  a  much  larger  wind 
shear  than  the  boundary  layer  over  flat  terrain.  The  total  integrated  shear  must  be  the 
same,  since  the  velocity  is  zero  at  the  wall  and  approaches  the  geostrophic  value  above 
the  PEL.  However,  the  shear  is  concentrated  in  a  shallow  layer  at  the  surface  over  flat 
terrain.  The  effect  on  particle  dispersion  is  therefore  limited  to  a  small  fraction  of  the 
cloud  near  the  ground.  The  profiles  over  the  ridges  show  that  the  shear  is  distributed  over 
a  much  deeper  layer,  certainly  extending  above  the  crests,  and  may  even  produce 
sigruficant  reverse  velocities  in  the  separation  regions.  This  increased  depth  of  mfluence 
is  responsible  for  the  increased  diffusion  in  the  x-direction,  since  a  large  fraction  of  the 
boundary  layer  is  now  affected  by  the  shear. 


We  have  not  been  able  to  perform  a  sufficient  number  of  dispersion  calculations 
to  provide  quantitative  statistics  on  cloud  spread  rates  for  the  range  of  meteorological 
conditions  and  release  locations  considered  in  the  study,  but  the  limited  information 
obtained  gives  insight  into  the  dominant  processes  that  need  to  considered.  A  simple 
means  of  estimating  the  dispersion  enhancement  over  complex  terrain  will  be  discussed 
in  Section  5.8. 
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5.7  PARTICLE  DEPOSITION  RATES. 

The  turbulent  deposition  of  small  particles  has  been  shown  to  be  an  important 
removal  mechanism  in  the  PBL,  and  practical  parameterization  schemes  were  discussed 
in  Section  4.  The  principal  flow  parameters  determining  the  deposition  rate  were  the 
surface  roughness  (or  vegetative  canopy  characteristics)  and  the  surface  stress  (or 
turbulence  levels).  The  presence  of  terrain  has  a  strong  influence  on  the  surface  stress 
and  local  turbulence  levels,  so  we  must  consider  the  effects  of  terrain  on  particle 

deposition  rates. 


The  T  FS  model  discussed  earlier  in  this  section  included  a  conservation  equation 
for  a  scalar  species  with  the  surface  deposition  parameterization  given  by  (4.18)-(4.21)  m 
Section  4.  Several  model  runs  were  conducted  with  a  scalar  species  introduced  into  the 
boundary  layer,  and  the  average  deposition  rate  over  the  computational  domain  was 
calculated.  The  rate  was  determined  as  the  ratio  of  the  average  surface  flux  of  scalar 
quantity  to  the  average  concentration  at  the  lowest  model  grid  level.  The  grid  location  is 
somewhat  arbitrary,  but  is  fixed  for  the  series  of  runs  and  therefore  provides  a  standard 
definition.  The  flux  of  scalar  depends  on  the  instantaneous,  local  value  of  the  surface 
stress,  as  does  the  local  concentration,  but  the  quantities  discussed  in  this  section  are 
averaged  over  the  4km  x  4km  domain  and  over  a  3000s  time  period. 

Table  5-2  shows  the  effective  deposition  velocities  for  a  number  of  LES  runs  with 
different  terrain  parameters.  It  is  apparent  from  the  table  that  the  slope  of  the  terrain  does 
not  produce  a  marked  effect  on  the  effective  deposition  velocity,  in  contrast  with  the 
momentum  transfer  results  of  Section  5.5.  It  seems  that  the  deposition  velocity  is  hardly 
changed  by  the  presence  of  terrain,  and  this  insensitivity  can  be  understood  from  the  fact 
that  the  turbulent  deposition  is  determined  by  the  surface  stress  distribution.  A  large 
component  of  the  increased  surface  force  is  due  to  the  pressure  force  on  the  hills,  while 
the  average  tangential  stress  component  is  generally  reduced  slightly  over  the  ridges.  The 
pressure  force  transfers  momentum  but  has  no  comparable  role  in  the  particle  deposition 
rate.  The  deposition  is  more  closely  related  to  the  tangential  force,  but  depends  on  the 
magnitude  of  the  stress.  The  jc-component  is  reduced  by  the  terrain ,  but  the  y-component 
is  generally  increased;  the  overall  average  magnitude  is  therefore  much  less  sensitive  to 

terrain  variations. 
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Table  5-2.  Deposition  velocity  for  lOpm  particles  over  terrain. 


Surface  roughness 
(m) 

Ridge  wavelength 
(km) 

Maximum  slope 

Deposition  velocity 
(cm/s) 

0.01 

2 

0.25 

0.4 

0.01 

2 

0.50 

0.6 

0.1 

2 

0.25 

1.4 

0.1 

2 

0.50 

1.5 

1.0 

1 

0.50 

5.7 

1.0 

2 

0.25 

4.8 

1.0 

2 

0.50 

6.3 

1.0 

4 

0.25 

6.7 

1.0 

4 

0.50 

6.3 

Table  5-2  strongly  suggest  that  the  representation  of  deposition  rates  over 
complex  terrain  should  utilize  the  tangential  force  parameterizations  from  (5.38)  in 
determining  the  effective  lu,  but  continue  to  use  the  flat  terrain  deposition  scheme  of 

Section  4. 


5.8  IMPLEMENTATION  IN  LATE-TIME  MODELS. 

We  have  briefly  suggested  some  of  the  directions  for  implementing  the  terrain 
representations  in  late-time  models.  One  of  the  key  aspects  of  such  representations  is  the 
need  to  describe  larger  scale  processes,  so  that  the  statistical  average  over  a  region  of 
complex  terrain  is  appropriate.  Detailed  calculations  for  particular  local  terrain 
conditions  requires  a  detailed  simulation  of  the  local  flow  field;  we  are  concerned  here 
with  larger  scale  effects.  We  need  to  represent  both  mean  boundary  layer  transport  and 
the  effective  dispersion  rates,  both  horizontal  and  vertical,  in  complex  terrain.  Deposition 
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rates  were  discussed  in  Section  5.7,  and  a  simple  scheme  was  recommended  based  on  the 
tangential  force  parameterization  (5.38). 

We  assume  that  the  large  scale  wind  field  used  to  drive  the  late-time  dispersion 
model  will  not  contain  a  detailed  terrain  parameterization  for  the  PBL,  so  that  the 
available  winds  are  closer  to  geostrophic  values.  Under  these  conditions  a  simple  time- 
dependent  scheme  for  estimating  the  mean  boundary  layer  wind  can  be  obtained  from  the 
integrated  momentum  equations  (5.25),  and  the  definition  of  the  mean  boundary  layer 
velocity,  mb,  in  (5.27).  The  equations  can  be  written 

^  =  -/(v.-v,)+^  (5-40a) 

or 

^  =  (S.-tOb) 

The  surface  force,  F,  is  determined  from  the  previously  developed 
parameterizations  (5.38),  and  the  inversion  depth,  Zi,  is  obtained  from  either  the 
meteorological  fields,  or  from  a  prediction  scheme  such  as  METPRO  (see  Section  2). 

The  dispersion  rates  in  complex  terrain  were  shown  to  be  a  result  of  both 
turbulence  and  mean  flow  effects.  The  turbulence  intensities  were  discussed  in  Section 
5.6,  where  the  surface  force  parameterization  was  suggested  as  a  basis  for  estimating  the 
large-scale  average  velocity  variance  profiles.  The  effect  of  the  mean  flow  perturbations 
can  be  modeled  by  including  a  mean  wind  shear  in  the  boundary  layer  Avind  profile.  The 
wind  shear  provides  a  persistent  effect  on  the  particle  dispersion,  and  models  the  effect  of 
the  slow  moving  air  in  the  sheltered  valley  regions.  We  emphasize  that  we  are 
representing  the  large  scale  processes  over  a  number  of  terrain  features,  not  the  detailed 
flow  distortion  around  a  particular  hill.  The  boundary  layer  wind  and  turbulence  profiles 
should  be  thought  of  as  a  large  scale  area  average,  over  a  horizontal  area  of  lOOkm^,  say. 

Unfortunately,  limitations  on  time  have  prevented  an  extensive  study  of  the 
quantitative  dispersion  mechanisms  for  the  range  of  meteorological  and  terrain 
parameters  studied  in  the  numerical  simulations.  We  therefore  suggest  a  basic 
parameterization  scheme,  with  numerical  parameters  specified  to  the  best  of  our  current 
knowledge.  We  expect  that  the  representation  could  be  refined  by  further  investigation  of 
the  LES  results.  The  assumption  in  the  scheme  is  that  the  effect  of  the  slow  moving  air  in 
the  valleys  can  be  modeled  as  a  simple  linear  shear  profile  in  the  velocity.  Under 
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homogeneous  mixing  conditions,  the  wind  shear  is  usually  confined  to  thin  layer  at  the 
ground  and  the  inversion.  The  shear  does  not  produce  strong  dispersion  of  material  under 
these  circumstances.  We  therefore  propose  that  the  wind  profile  include  a  linear  regime 
in  the  lower  part  of  the  boundary  layer,  with  a  maximum  reduction  that  depends  on  the 
r.m.s.  slope  of  the  topography,  and  a  depth  that  depends  on  the  r.m.s.  height  of  the  terrain. 
The  simplest  model  incorporating  these  effects  is 


Ui=UQi-C2 


1-^ 


for  z<ciOh,  and  u=uo  above.  The  mean  velocity  is  required  to  be  mb.  so  this  determines  uq 
if  the  other  quantities  are  prescribed.  As  a  preliminary  estimate,  we  suggest  that  ci=2.5, 

and  C2=2. 


SECTION  6 
CONCLUSIONS 

The  representation  of  planetary  boundary  effects  on  the  distribution  of  near¬ 
surface  dust  in  late-time  modeling  of  nuclear  clouds  has  been  considered. 
Representations  of  the  boundary  layer  turbulence  and  effective  diffusivity  have  been 
recommended  for  application  in  the  models.  The  diffusivity  estimate  is  based  on  recent 
turbulence  closure  models,  and  provides  a  generalized  description  of  the  boundary  layer 
based  on  a  small  number  of  boundary  layer  parameters.  The  boundary  layer  parameters 
include  the  wind  speed  and  surface  fluxes,  in  addition  to  the  surface  roughness  and 
boundary  layer  depth,  and  a  simplified  scheme  for  determining  the  fluxes  using  the 
METPRO  meteorological  pre-processor  (Paine,  1987)  is  recommended. 

A  relatively  detailed  examination  of  the  effect  of  dust  on  the  solar  radiative  flux 
transport  was  conducted,  since  the  solar  heating  is  a  major  component  of  the  boundary 
layer  energy  input  A  preliminary  study  indicated  that  concentrations  of  small  particles, 
i.e.,  smaller  than  10pm  diameter,  at  levels  of  around  10"‘^g/cc  could  produce  significant 
effects  on  the  radiative  transfer.  Such  concentrations  are  feasible  in  a  multiburst  scenario, 
so  detailed  dynamic  calculations  were  performed  using  the  ARAP  turbulence  closure 
model  and  the  delta-Eddington  radiative  flux  model.  The  effect  of  a  dense  dust  cloud  is 
to  absorb  the  solar  flux  in  the  upper  portion  of  the  cloud,  preventing  the  flux  from 
reaching  the  ground.  This  removes  the  buoyancy  generation  mechanism  from  the 
planetary  boundary  layer,  and  suppresses  the  turbulent  diffusion  near  the  surface.  The 
lower  part  of  the  cloud  remains  unmixed  during  the  daytime  as  a  result  of  this 
suppression.  The  absorbing  part  of  the  dust  cloud,  however,  is  lofted  by  the  heating  and 
maintains  a  low  level  of  turbulent  mixing.  The  representation  of  these  effects  in  a  late¬ 
time  dispersion  model  is  a  difficult  problem.  Late-time  models  are  traditionally  passive 
in  nature,  i.e.,  the  dust  is  transported  and  diffused  but  does  not  induce  its  own  velocity 
field.  The  effective  radiation  reaching  the  surface  can  be  estimated  from  area-averaged 
vertical  integrals  of  the  dust,  so  that  a  crude  estimate  of  the  surface  heating  changes  can 
be  included,  but  the  representation  of  the  dynamic  lofting  mechanism  requires  a  more 
detailed  dynamic  model 
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Turbulent  deposition  of  particles  at  the  surface  is  an  important  process  for  small 
dust  particles,  and  parameterizations  for  inclusion  in  late-time  models  have  been 
suggested.  The  models  represent  the  effects  of  particle  size  and  the  nature  of  the 
vegetative  canopy  or  surface  roughness,  if  details  are  available.  The  models  are  simple 
enough  for  practical  calculation,  and  provide  a  reasonably  accurate  description  of 
laboratory  data. 

Most  of  the  existing  boundary  layer  representations  are  appropriate  for  flat, 
homogeneous  terrain,  but  most  real  applications  involve  flow  over  hills  and  valleys.  In 
order  to  improve  our  understanding  of  flow  over  terrain  and  provide  a  basis  for 
representing  the  effects  in  late-time  dispersion  models,  numerical  simulations  of 
atmospheric  boundary  layer  flow  over  topography  under  convective  conditions  have  been 
performed  for  a  range  of  meteorological  conditions  and  terrain  shapes.  The  bulk  of  the 
turbulent  velocity  fluctuations  have  been  resolved  using  the  Large-Eddy  Simulation 
technique  with  a  terrain-following  coordinate  transformation.  An  extensive  set  of 
simulations  was  conducted  to  examine  the  flow  response  to  variations  in  terrain 
amplitude  and  wavelength,  surface  roughness,  boundary  layer  stability,  and  geostrophic 
flow  direction  for  idealized  sinusoidal  ridges. 

Terrain  wavelengths  between  1km  and  4km  were  studied,  with  maximum  slopes 
up  to  0.5.  Surface  roughness  length  was  varied  between  1cm  and  Im,  and  geostrophic 
flow  direction  was  varied  relative  to  the  ridge  orientation.  In  addition,  the  stability  was 
varied  by  changing  the  surface  heat  flux  and  geostrophic  wind  speed.  Mean  flow  Eelds 
and  surface  forces  were  the  principal  focus  of  the  study.  The  surface  forces  were 
analyzed  and  used  to  develop  a  simplified  parameterization  scheme  for  late-time 
dispersion  models  using  the  average  boundary  layer  velocity. 

The  LES  results  have  been  used  to  derive  a  simplified  surface  drag 
parameterization  for  convective  flow  over  more  general  topography,  utilizing  the  results 
from  simulations  with  different  shaped  ridges  and  three-dimensionality.  The  simulations 
also  provide  insight  into  the  enhanced  horizontal  dispersion  rates  in  complex  terrain,  and 
a  simple  estimate  for  general  conditions  has  been  recommended.  The  LES  calculations 
showed  that  particle  deposition  rates  are  not  as  strongly  modified  by  terrain,  and  the  flat 
surface  parameterization  can  be  used  in  complex  terrain,  provided  the  surface  friction 
velocity  is  estimated  from  the  tangential  stress  force,  i.e.,  the  pressure  force  should  not  be 
included. 
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